## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.2 ──
## ✔ ggplot2 3.3.6 ✔ purrr 0.3.4
## ✔ tibble 3.2.1 ✔ dplyr 1.1.2
## ✔ tidyr 1.2.0 ✔ stringr 1.4.1
## ✔ readr 2.1.2 ✔ forcats 0.5.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
##
## Attaching package: 'lubridate'
##
##
## The following objects are masked from 'package:base':
##
## date, intersect, setdiff, union
##
##
## corrplot 0.92 loaded
##
##
## Attaching package: 'ggpp'
##
##
## The following object is masked from 'package:ggplot2':
##
## annotate
##
##
##
## Attaching package: 'flextable'
##
##
## The following object is masked from 'package:purrr':
##
## compose
##
##
##
## Attaching package: 'maps'
##
##
## The following object is masked from 'package:purrr':
##
## map
##
##
## ℹ Google's Terms of Service: ]8;;https://mapsplatform.google.com<https://mapsplatform.google.com>]8;;
## ℹ Please cite ggmap if you use it! Use `citation("ggmap")` for details.
##
## Attaching package: 'scales'
##
##
## The following object is masked from 'package:purrr':
##
## discard
##
##
## The following object is masked from 'package:readr':
##
## col_factor
##
##
## Loading required package: sp
##
## Checking rgeos availability: TRUE
## Please note that 'maptools' will be retired by the end of 2023,
## plan transition at your earliest convenience;
## some functionality will be moved to 'sp'.
##
##
## Attaching package: 'raster'
##
##
## The following object is masked from 'package:dplyr':
##
## select
##
##
## rgeos version: 0.6-3, (SVN revision 696)
## GEOS runtime version: 3.10.2-CAPI-1.16.0
## Please note that rgeos will be retired during October 2023,
## plan transition to sf or terra functions using GEOS at your earliest convenience.
## See https://r-spatial.org/r/2023/05/15/evolution4.html for details.
## GEOS using OverlayNG
## Linking to sp version: 1.6-0
## Polygon checking: TRUE
##
##
##
## Attaching package: 'rgeos'
##
##
## The following object is masked from 'package:dplyr':
##
## symdiff
This script reads in cleaned datafiles with stream temperature and invertebrate data produced by the invert_data_cleaning.Rmd and teton_data_cleaning.Rmd scripts and performs more involved analyses of stream temperature and invertebrate density/diversity over time. A few findings of interest from previous analyses are that:
Average summer (July-Aug) stream temperatures have not changed significantly over time
Average summer (July-Aug) air temperatures have increased over time
Overall, invetebrate richness has not changed over time
Invertebrate communities in snow-fed streams appear to be more variable than those in glacier-fed streams
To build on these findings, the main sections of this script are:
I. Temperature:
Trends in summer mean, max, min, and degree-days over time: Site-by-site and grouped by source
Modeling stream temperature using SNOTEL data: Site-by-site and grouped by source
Linear mixed models of temperature over time
II. Invertebrates:
Trends in richness, abundance, and density over time: Site-by-site and grouped by source
Modeling richness, abundance, and density using SNOTEL data: Site-by-site and grouped by source
III. Temperature and Invertebrates
Read in necessary files:
stream_temp <- read.csv("Temperature/cleaned_full_datasets/temps_hourly.csv") #hourly stream temp data
invert_density <- read.csv("invert_data/cleaned_csv/full_invert_densities.csv")%>%
rename(site = Stream)#invert density data
invert_richness <- read.csv("invert_data/cleaned_csv/invert_richness.csv") #invert richness data
snotel <- read.csv("teton_snotel.csv") #snotel data
sources <- read.csv("source_info.csv")%>%
rename(site = stream) #source info, rename for merge
### adding source info to stream temps and invert datasets:
stream_source <- merge(stream_temp, sources, all = T)%>%
mutate(date1 = ymd(date1),
year = year(date1))%>%
filter(!is.na(temp_c))#Add source info to hourly temperature data; re-code date as r date format
stream_source_daily <- stream_source%>%
group_by(site, full_name, stream_code, date1, year, lat, long)%>%
summarise(t_xbar = mean(temp_c, na.rm = T),
t_max = max(temp_c, na.rm = T),
t_min = min(temp_c, na.rm = T),
t_range = t_max-t_min,
source = unique(source)) #calculate daily mean, min, and max temps for each site along with temp range
## `summarise()` has grouped output by 'site', 'full_name', 'stream_code',
## 'date1', 'year', 'lat'. You can override using the `.groups` argument.
density_source <- merge(invert_density, sources, all = T)
richness_source <- merge(invert_richness, sources, all = T)
First, calculate summer (july-aug) mean daily temp, avg daily max and min; plus degree days over 15c? for each site:
stream_summer <- stream_source_daily %>%
mutate(mo = month(date1))%>% #add month column
filter(mo == 7 | mo == 8)%>% #extract July and August
group_by(site, full_name, stream_code, year, lat, long)%>% #group by site and year
summarise(summer_xbar = mean(t_xbar, na.rm = T),
max_xbar = mean(t_max, na.rm = T),
min_xbar = mean(t_min, na.rm= T), #calculate mean temperature, mean daily max and min
range_xbar = mean(t_range, na.rm = T), #calculate mean daily temperature range
dd_2.5 = length(which(t_max >= 2.5)),
dd_5 = length(which(t_max >= 5)),
dd_10 = length(which(t_max >= 10)), #calculate no. days with max temp >2.5, 5, and 10 deg
source = unique(source)) #retain source col
## `summarise()` has grouped output by 'site', 'full_name', 'stream_code', 'year',
## 'lat'. You can override using the `.groups` argument.
Next, extract sites with at least 3 years of data:
under3 <- stream_summer %>%
group_by(site)%>% #group by site
filter(!is.na(year))%>% #remove rows with no year
summarise(n_yrs = length(unique(year)))%>% #count num unique years for site
filter(n_yrs < 3) #extract sites with <3yrs data
summer_3plus <- stream_summer %>%
filter(!site%in%under3$site)
Trends in mean temp:
Nest data by site:
summer_nested <- summer_3plus %>%
filter(!is.na(year))%>%
nest(data = -site)
lm for each site:
summer.mean.lms <- summer_nested %>%
mutate(model = purrr::map(data, ~lm(summer_xbar~year, data = .)%>% #run a model of temp~year for each site (referenced by "."); store results in new column called "model"
tidy))%>% #tidy function to construct a summary table with model results (seperate tables for each site)
unnest(model) %>% #unnest model column to get rows for each site, columns for model metrics
filter(term == "year") #remove intercept term to just get info on year for each site
summer.mean.lms
## # A tibble: 12 × 7
## # Groups: site [12]
## site data term estimate std.error statistic p.value
## <chr> <list> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 cloudveil <gropd_df [4 × 13]> year 0.0500 0.265 0.188 0.868
## 2 delta <gropd_df [6 × 13]> year 0.0321 0.0321 0.999 0.374
## 3 grizzly <gropd_df [5 × 13]> year 0.198 0.561 0.353 0.747
## 4 mid_teton <gropd_df [8 × 13]> year 0.00372 0.0258 0.144 0.890
## 5 n_teton <gropd_df [8 × 13]> year -0.137 0.308 -0.444 0.673
## 6 paintbrush <gropd_df [6 × 13]> year 0.442 0.471 0.937 0.402
## 7 s_ak_basin <gropd_df [6 × 13]> year -0.0169 0.0966 -0.175 0.870
## 8 s_cascade <gropd_df [5 × 13]> year -0.151 0.147 -1.03 0.379
## 9 s_teton <gropd_df [6 × 13]> year 0.153 0.354 0.433 0.688
## 10 silver_run <gropd_df [5 × 13]> year 0.419 0.0856 4.90 0.0163
## 11 skillet <gropd_df [5 × 13]> year 0.165 0.205 0.807 0.479
## 12 windcave <gropd_df [5 × 13]> year 0.102 0.0502 2.03 0.135
Nothing significant at 0.05 level; weakly significant increase in mean summer temp at S AK basin (p = 0.08)
Visualization:
summer_mean_withps <- merge(summer_3plus, summer.mean.lms)
pres.pal <- c("#C26ED6", "#F89225", "#76D96F")
summer.mean <- ggplot(summer_mean_withps, aes(x = year, y = summer_xbar, color = source))+
geom_point()+
geom_line()+
geom_smooth(se = F, method = "lm", linetype = "dashed", size = 0.5)+
theme_classic()+
facet_wrap(~full_name, scales = "free")+
scale_color_manual(values = pres.pal, labels = c("Glacier", "Snowmelt", "Icy Seep"))+
geom_text_npc(aes(npcx = 0.9, npcy = 0.1, label = paste( "P = ",round(p.value, 2), sep = "")))+
labs(x = "Date", y = "Mean daily July-Aug. stream temp., C", color = "Stream Source")+
theme(legend.position = "bottom",
axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
axis.title = element_text(size = 11, face = "bold"),
axis.text = element_text(size = 8),
legend.text = element_text(size = 10),
legend.title = element_text(size = 11, face = "bold"),
strip.text = element_text(size = 10, face = "bold"))
summer.mean
## `geom_smooth()` using formula 'y ~ x'
Mean daily max temp
summer.max.lms <- summer_nested %>%
mutate(model = purrr::map(data, ~lm(max_xbar~year, data = .)%>% #run a model of temp~year for each site (referenced by "."); store results in new column called "model"
tidy))%>% #tidy function to construct a summary table with model results (seperate tables for each site)
unnest(model) %>% #unnest model column to get rows for each site, columns for model metrics
filter(term == "year") #remove intercept term to just get info on year for each site
summer.max.lms
## # A tibble: 12 × 7
## # Groups: site [12]
## site data term estimate std.error statistic p.value
## <chr> <list> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 cloudveil <gropd_df [4 × 13]> year 0.0266 0.314 0.0844 0.940
## 2 delta <gropd_df [6 × 13]> year 0.0781 0.107 0.729 0.506
## 3 grizzly <gropd_df [5 × 13]> year -0.0597 0.915 -0.0653 0.952
## 4 mid_teton <gropd_df [8 × 13]> year -0.0808 0.0565 -1.43 0.203
## 5 n_teton <gropd_df [8 × 13]> year -0.320 0.395 -0.810 0.449
## 6 paintbrush <gropd_df [6 × 13]> year 0.694 0.774 0.897 0.420
## 7 s_ak_basin <gropd_df [6 × 13]> year -0.00900 0.185 -0.0487 0.963
## 8 s_cascade <gropd_df [5 × 13]> year -0.353 0.403 -0.877 0.445
## 9 s_teton <gropd_df [6 × 13]> year 0.00464 0.552 0.00841 0.994
## 10 silver_run <gropd_df [5 × 13]> year 0.293 0.0926 3.16 0.0509
## 11 skillet <gropd_df [5 × 13]> year -0.0680 0.380 -0.179 0.869
## 12 windcave <gropd_df [5 × 13]> year 0.140 0.0668 2.10 0.127
Significant increase at Skillet, non-significant increases at all but S cascade and mid teton
Visualization:
summer_max_withps <- merge(summer_3plus, summer.max.lms)
summer.max <- ggplot(summer_max_withps, aes(x = year, y = max_xbar, color = source))+
geom_point()+
geom_line()+
geom_smooth(se = F, method = "lm", linetype = "dashed", size = 0.5)+
theme_classic()+
facet_wrap(~full_name, scales = "free")+
scale_color_manual(values = pres.pal, labels = c("Glacier", "Snowmelt", "Icy Seep"))+
geom_text_npc(aes(npcx = 0.9, npcy = 0.1, label = paste( "P = ",round(p.value, 2), sep = "")))+
labs(x = "Date", y = "Mean daily max. July-Aug. stream temp., C", color = "Stream Source")+
theme(legend.position = "bottom",
axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
axis.title = element_text(size = 11, face = "bold"),
axis.text = element_text(size = 8),
legend.text = element_text(size = 10),
legend.title = element_text(size = 11, face = "bold"),
strip.text = element_text(size = 10, face = "bold"))
summer.max
## `geom_smooth()` using formula 'y ~ x'
Mean daily min temp
summer.min.lms <- summer_nested %>%
mutate(model = purrr::map(data, ~lm(min_xbar~year, data = .)%>% #run a model of temp~year for each site (referenced by "."); store results in new column called "model"
tidy))%>% #tidy function to construct a summary table with model results (seperate tables for each site)
unnest(model) %>% #unnest model column to get rows for each site, columns for model metrics
filter(term == "year") #remove intercept term to just get info on year for each site
summer.min.lms
## # A tibble: 12 × 7
## # Groups: site [12]
## site data term estimate std.error statistic p.value
## <chr> <list> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 cloudveil <gropd_df [4 × 13]> year 0.0663 0.223 0.297 0.794
## 2 delta <gropd_df [6 × 13]> year 0.0159 0.00768 2.07 0.107
## 3 grizzly <gropd_df [5 × 13]> year 0.252 0.363 0.695 0.537
## 4 mid_teton <gropd_df [8 × 13]> year 0.0480 0.0186 2.59 0.0414
## 5 n_teton <gropd_df [8 × 13]> year -0.0138 0.227 -0.0607 0.954
## 6 paintbrush <gropd_df [6 × 13]> year 0.284 0.307 0.926 0.407
## 7 s_ak_basin <gropd_df [6 × 13]> year -0.0203 0.0714 -0.285 0.790
## 8 s_cascade <gropd_df [5 × 13]> year -0.0329 0.0169 -1.95 0.146
## 9 s_teton <gropd_df [6 × 13]> year 0.206 0.215 0.955 0.394
## 10 silver_run <gropd_df [5 × 13]> year 0.433 0.139 3.13 0.0522
## 11 skillet <gropd_df [5 × 13]> year 0.268 0.107 2.51 0.0872
## 12 windcave <gropd_df [5 × 13]> year 0.0867 0.0350 2.48 0.0892
Non-significant increases at all sites except south cascade
Visualization:
summer_min_withps <- merge(summer_3plus, summer.min.lms)
summer.min <- ggplot(summer_min_withps, aes(x = year, y = min_xbar, color = source))+
geom_point()+
geom_line()+
geom_smooth(se = F, method = "lm", linetype = "dashed", size = 0.5)+
theme_classic()+
facet_wrap(~full_name, scales = "free")+
scale_color_manual(values = pres.pal, labels = c("Glacier", "Snowmelt", "Icy Seep"))+
geom_text_npc(aes(npcx = 0.9, npcy = 0.1, label = paste( "P = ",round(p.value, 2), sep = "")))+
labs(x = "Date", y = "Mean daily min. July-Aug. stream temp., C", color = "Stream Source")+
theme(legend.position = "bottom",
axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
axis.title = element_text(size = 11, face = "bold"),
axis.text = element_text(size = 8),
legend.text = element_text(size = 10),
legend.title = element_text(size = 11, face = "bold"),
strip.text = element_text(size = 10, face = "bold"))
summer.min
## `geom_smooth()` using formula 'y ~ x'
Daily temperature fluctuation
summer.range.lms <- summer_nested %>%
mutate(model = purrr::map(data, ~lm(range_xbar~year, data = .)%>% #run a model of temp~year for each site (referenced by "."); store results in new column called "model"
tidy))%>% #tidy function to construct a summary table with model results (seperate tables for each site)
unnest(model) %>% #unnest model column to get rows for each site, columns for model metrics
filter(term == "year") #remove intercept term to just get info on year for each site
summer.range.lms
## # A tibble: 12 × 7
## # Groups: site [12]
## site data term estimate std.error statistic p.value
## <chr> <list> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 cloudveil <gropd_df [4 × 13]> year -0.0398 0.0978 -0.407 0.724
## 2 delta <gropd_df [6 × 13]> year 0.0622 0.109 0.573 0.597
## 3 grizzly <gropd_df [5 × 13]> year -0.312 0.573 -0.545 0.624
## 4 mid_teton <gropd_df [8 × 13]> year -0.129 0.0493 -2.61 0.0399
## 5 n_teton <gropd_df [8 × 13]> year -0.306 0.171 -1.79 0.123
## 6 paintbrush <gropd_df [6 × 13]> year 0.410 0.496 0.826 0.455
## 7 s_ak_basin <gropd_df [6 × 13]> year 0.0113 0.172 0.0659 0.951
## 8 s_cascade <gropd_df [5 × 13]> year -0.321 0.389 -0.823 0.471
## 9 s_teton <gropd_df [6 × 13]> year -0.201 0.342 -0.587 0.588
## 10 silver_run <gropd_df [5 × 13]> year -0.141 0.203 -0.695 0.537
## 11 skillet <gropd_df [5 × 13]> year -0.336 0.276 -1.22 0.310
## 12 windcave <gropd_df [5 × 13]> year 0.0535 0.0346 1.55 0.220
Non-significant increases at all sites except south cascade
Visualization:
summer_range_withps <- merge(summer_3plus, summer.range.lms)
summer.min <- ggplot(summer_range_withps, aes(x = year, y = range_xbar, color = source))+
geom_point()+
geom_line()+
geom_smooth(se = F, method = "lm", linetype = "dashed", size = 0.5)+
theme_classic()+
facet_wrap(~full_name, scales = "free")+
scale_color_manual(values = pres.pal, labels = c("Glacier", "Snowmelt", "Icy Seep"))+
geom_text_npc(aes(npcx = 0.9, npcy = 0.1, label = paste( "P = ",round(p.value, 2), sep = "")))+
labs(x = "Date", y = "Mean daily July-Aug. stream temp. change, C", color = "Stream Source")+
theme(legend.position = "bottom",
axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
axis.title = element_text(size = 11, face = "bold"),
axis.text = element_text(size = 8),
legend.text = element_text(size = 10),
legend.title = element_text(size = 11, face = "bold"),
strip.text = element_text(size = 10, face = "bold"))
summer.min
## `geom_smooth()` using formula 'y ~ x'
Degree Days
Days over 2.5 C:
summer.2.5.lms <- summer_nested %>%
mutate(model = purrr::map(data, ~lm(dd_2.5~year, data = .)%>% #run a model of temp~year for each site (referenced by "."); store results in new column called "model"
tidy))%>% #tidy function to construct a summary table with model results (seperate tables for each site)
unnest(model) %>% #unnest model column to get rows for each site, columns for model metrics
filter(term == "year") #remove intercept term to just get info on year for each site
summer.2.5.lms
## # A tibble: 12 × 7
## # Groups: site [12]
## site data term estimate std.error statistic p.value
## <chr> <list> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 cloudveil <gropd_df [4 × 13]> year 2.00 6.81 0.294 0.797
## 2 delta <gropd_df [6 × 13]> year 2.89 3.74 0.772 0.483
## 3 grizzly <gropd_df [5 × 13]> year 2.50 4.14 0.604 0.588
## 4 mid_teton <gropd_df [8 × 13]> year 2.94 2.41 1.22 0.269
## 5 n_teton <gropd_df [8 × 13]> year 3.18 1.78 1.78 0.125
## 6 paintbrush <gropd_df [6 × 13]> year 4.34 2.21 1.97 0.121
## 7 s_ak_basin <gropd_df [6 × 13]> year 0.149 1.34 0.111 0.917
## 8 s_cascade <gropd_df [5 × 13]> year 0.407 2.33 0.175 0.872
## 9 s_teton <gropd_df [6 × 13]> year 3.09 3.73 0.828 0.454
## 10 silver_run <gropd_df [5 × 13]> year -1.60 5.42 -0.295 0.787
## 11 skillet <gropd_df [5 × 13]> year 0.700 4.26 0.164 0.880
## 12 windcave <gropd_df [5 × 13]> year 0.948 1.47 0.646 0.564
Non-significant increases at all sites except south cascade
Visualization:
summer_2.5_withps <- merge(summer_3plus, summer.2.5.lms)
summer.2.5 <- ggplot(summer_2.5_withps, aes(x = year, y = dd_2.5, color = source))+
geom_point()+
geom_line()+
geom_smooth(se = F, method = "lm", linetype = "dashed", size = 0.5)+
theme_classic()+
facet_wrap(~full_name, scales = "free")+
scale_color_manual(values = pres.pal, labels = c("Glacier", "Snowmelt", "Icy Seep"))+
geom_text_npc(aes(npcx = 0.9, npcy = 0.1, label = paste( "P = ",round(p.value, 2), sep = "")))+
labs(x = "Date", y = "No. summer (July-Aug.) days with max temp >2.5 C", color = "Stream Source")+
theme(legend.position = "bottom",
axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
axis.title = element_text(size = 11, face = "bold"),
axis.text = element_text(size = 8),
legend.text = element_text(size = 10),
legend.title = element_text(size = 11, face = "bold"),
strip.text = element_text(size = 10, face = "bold"))
summer.2.5
## `geom_smooth()` using formula 'y ~ x'
Days over 5 C:
summer.5.lms <- summer_nested %>%
mutate(model = purrr::map(data, ~lm(dd_5~year, data = .)%>% #run a model of temp~year for each site (referenced by "."); store results in new column called "model"
tidy))%>% #tidy function to construct a summary table with model results (seperate tables for each site)
unnest(model) %>% #unnest model column to get rows for each site, columns for model metrics
filter(term == "year") #remove intercept term to just get info on year for each site
summer.5.lms
## # A tibble: 12 × 7
## # Groups: site [12]
## site data term estimate std.error statistic p.value
## <chr> <list> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 cloudveil <gropd_df [4 × 13]> year 0.5 1.32 0.378 0.742
## 2 delta <gropd_df [6 × 13]> year 0 0 NaN NaN
## 3 grizzly <gropd_df [5 × 13]> year 2.30 4.05 0.568 0.610
## 4 mid_teton <gropd_df [8 × 13]> year -0.0238 0.117 -0.203 0.846
## 5 n_teton <gropd_df [8 × 13]> year 2.55 1.63 1.56 0.169
## 6 paintbrush <gropd_df [6 × 13]> year 5.97 1.81 3.30 0.0299
## 7 s_ak_basin <gropd_df [6 × 13]> year 0.106 0.751 0.141 0.895
## 8 s_cascade <gropd_df [5 × 13]> year -0.923 2.01 -0.458 0.678
## 9 s_teton <gropd_df [6 × 13]> year 2.97 3.67 0.810 0.464
## 10 silver_run <gropd_df [5 × 13]> year -1.60 5.29 -0.303 0.782
## 11 skillet <gropd_df [5 × 13]> year 0.100 3.42 0.0292 0.979
## 12 windcave <gropd_df [5 × 13]> year 0.955 0.577 1.66 0.196
Non-significant increases at all sites except south cascade, mid teton, and delta (0 all years)
Visualization:
summer_5_withps <- merge(summer_3plus, summer.5.lms)
summer.5 <- ggplot(summer_5_withps, aes(x = year, y = dd_5, color = source))+
geom_point()+
geom_line()+
geom_smooth(se = F, method = "lm", linetype = "dashed", size = 0.5)+
theme_classic()+
facet_wrap(~full_name, scales = "free")+
scale_color_manual(values = pres.pal, labels = c("Glacier", "Snowmelt", "Icy Seep"))+
geom_text_npc(aes(npcx = 0.9, npcy = 0.1, label = paste( "P = ",round(p.value, 2), sep = "")))+
labs(x = "Date", y = "No. summer (July-Aug.) days with max temp >5 C", color = "Stream Source")+
theme(legend.position = "bottom",
axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
axis.title = element_text(size = 11, face = "bold"),
axis.text = element_text(size = 8),
legend.text = element_text(size = 10),
legend.title = element_text(size = 11, face = "bold"),
strip.text = element_text(size = 10, face = "bold"))
summer.5
## `geom_smooth()` using formula 'y ~ x'
Days over 10 C:
summer.10.lms <- summer_nested %>%
mutate(model = purrr::map(data, ~lm(dd_10~year, data = .)%>% #run a model of temp~year for each site (referenced by "."); store results in new column called "model"
tidy))%>% #tidy function to construct a summary table with model results (seperate tables for each site)
unnest(model) %>% #unnest model column to get rows for each site, columns for model metrics
filter(term == "year") #remove intercept term to just get info on year for each site
summer.10.lms
## # A tibble: 12 × 7
## # Groups: site [12]
## site data term estimate std.error statistic p.value
## <chr> <list> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 cloudveil <gropd_df [4 × 13]> year 0 0 NaN NaN
## 2 delta <gropd_df [6 × 13]> year 0 0 NaN NaN
## 3 grizzly <gropd_df [5 × 13]> year 1.40 4.11 0.341 0.756
## 4 mid_teton <gropd_df [8 × 13]> year 0 0 NaN NaN
## 5 n_teton <gropd_df [8 × 13]> year 1.82 1.72 1.06 0.332
## 6 paintbrush <gropd_df [6 × 13]> year 2.71 1.28 2.13 0.100
## 7 s_ak_basin <gropd_df [6 × 13]> year -0.0124 0.176 -0.0705 0.947
## 8 s_cascade <gropd_df [5 × 13]> year 0 0 NaN NaN
## 9 s_teton <gropd_df [6 × 13]> year 2.34 3.55 0.661 0.545
## 10 silver_run <gropd_df [5 × 13]> year 0.700 4.70 0.149 0.891
## 11 skillet <gropd_df [5 × 13]> year 0 0 NaN NaN
## 12 windcave <gropd_df [5 × 13]> year 0 0 NaN NaN
Nothing significant.
Visualization:
summer_10_withps <- merge(summer_3plus, summer.10.lms)
summer.10 <- ggplot(summer_10_withps, aes(x = year, y = dd_10, color = source))+
geom_point()+
geom_line()+
geom_smooth(se = F, method = "lm", linetype = "dashed", size = 0.5)+
theme_classic()+
facet_wrap(~full_name, scales = "free")+
scale_color_manual(values = pres.pal, labels = c("Glacier", "Snowmelt", "Icy Seep"))+
geom_text_npc(aes(npcx = 0.9, npcy = 0.1, label = paste( "P = ",round(p.value, 2), sep = "")))+
labs(x = "Date", y = "No. summer (July-Aug.) days with max temp >10 C", color = "Stream Source")+
theme(legend.position = "bottom",
axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
axis.title = element_text(size = 11, face = "bold"),
axis.text = element_text(size = 8),
legend.text = element_text(size = 10),
legend.title = element_text(size = 11, face = "bold"),
strip.text = element_text(size = 10, face = "bold"))
summer.10
## `geom_smooth()` using formula 'y ~ x'
Make and save a table with stream names, years of data available for each, and stream source:
First, table set-up:
stream_info <- as.data.frame(summer_3plus) %>%
dplyr::select(full_name, year, source)%>% #dplyr::select full name, year, and source columns
group_by(full_name, source)%>% #group by full name and source
summarise(range = paste(min(year, na.rm = T), max(year, na.rm = T), sep = "-"))%>% #add date range column for each stream
arrange(source, full_name)%>% #sort by source then stream name
mutate(source = case_when(source == "glacier"~"Glacier", source == "snowmelt"~"Snow melt", source == "sub_ice"~"Subterranean ice"))
## `summarise()` has grouped output by 'full_name'. You can override using the
## `.groups` argument.
Convert to FlexTable and export:
tbl1 <- flextable(stream_info)%>%
theme_vanilla()%>%
align(part = "all", align = "left")%>%
set_header_labels(full_name = "Stream Name", source = "Hydrologic Source", range = "Date Range")%>%
set_table_properties(layout = "autofit")
tbl1
Stream Name | Hydrologic Source | Date Range |
|---|---|---|
Cloudveil | Glacier | 2019-2022 |
Delta | Glacier | 2017-2022 |
Middle Teton | Glacier | 2015-2022 |
Skillet | Glacier | 2018-2022 |
Grizzly | Snow melt | 2018-2022 |
North Teton | Snow melt | 2015-2022 |
Paintbrush | Snow melt | 2017-2022 |
Silver Run | Snow melt | 2017-2021 |
South Teton | Snow melt | 2017-2022 |
S. AK Basin | Subterranean ice | 2016-2022 |
South Cascade | Subterranean ice | 2015-2022 |
Wind Cave | Subterranean ice | 2015-2021 |
save_as_docx(tbl1,path = "manuscript/tables/table1_streaminfo.docx", align = "center")
save_as_image(tbl1,path = "manuscript/tables/table1_streaminfo.png", res = 300)
## [1] "manuscript/tables/table1_streaminfo.png"
#Data set-up:
mean.long <- mean(summer_3plus$long, na.rm=T)
mean.lat <- mean(summer_3plus$lat, na.rm =T) #calculate mean lat/long of all sites
map_dat <- as.data.frame(summer_3plus)%>%
dplyr::select(full_name, source, lat, long)%>%
distinct()%>%
filter(!is.na(lat))
# Basemap: use get_map to pull map from google API:
basemap <- get_map(c(mean.long,mean.lat), #center map on mean lat and long of all sites
source = "stamen", zoom = 11, maptype = "terrain-background") #set map type and zoom
## ℹ <]8;;https://maps.googleapis.com/maps/api/staticmap?center=43.742633,-110.825418&zoom=11&size=640x640&scale=2&maptype=terrain&key=xxxhttps://maps.googleapis.com/maps/api/staticmap?center=43.742633,-110.825418&zoom=11&size=640x640&scale=2&maptype=terrain&key=xxx]8;;>
## ℹ Map tiles by Stamen Design, under CC BY 3.0. Data by OpenStreetMap, under ODbL.
bb <- attr(basemap, "bb") #extract max/min lat and long of the basemap for clipping
##GTNP boundary:
gtnp <- rgdal::readOGR("gis_data/gnp_gtnp")%>% #read in stored shapefile of GTNP and GNP polygons
spTransform(crs("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")) #Update CRS to lat/long (comes in as UTM)
## Please note that rgdal will be retired during October 2023,
## plan transition to sf/stars/terra functions using GDAL and PROJ
## at your earliest convenience.
## See https://r-spatial.org/r/2023/05/15/evolution4.html and https://github.com/r-spatial/evolution
## rgdal: version: 1.6-7, (SVN revision 1203)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.4.2, released 2022/03/08
## Path to GDAL shared files: /Users/GordonG/Library/R/x86_64/4.2/library/rgdal/gdal
## GDAL binary built with GEOS: FALSE
## Loaded PROJ runtime: Rel. 8.2.1, January 1st, 2022, [PJ_VERSION: 821]
## Path to PROJ shared files: /Users/GordonG/Library/R/x86_64/4.2/library/rgdal/proj
## PROJ CDN enabled: FALSE
## Linking to sp version:1.6-1
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading sp or rgdal.
## Warning: OGR support is provided by the sf and terra packages among others
## Warning: OGR support is provided by the sf and terra packages among others
## Warning: OGR support is provided by the sf and terra packages among others
## Warning: OGR support is provided by the sf and terra packages among others
## Warning: OGR support is provided by the sf and terra packages among others
## Warning: OGR support is provided by the sf and terra packages among others
## Warning: OGR support is provided by the sf and terra packages among others
## OGR data source with driver: ESRI Shapefile
## Source: "/Users/GordonG/iCloud Drive (Archive)/Documents/Work/hotaling_data_analyist/tasr/teton_trends/gis_data/gnp_gtnp", layer: "gnp_gtnp"
## with 2 features
## It has 20 fields
## Warning: PROJ support is provided by the sf and terra packages among others
gtnp_clip <- crop(gtnp, extent(bb[1,2], bb[1,4], bb[1,1], bb[1,3]))%>% #Clip GTNP shapefile to the extent of the basemap
fortify() #fortify to convert from SpatialPolygon to dataframe for plotting
## Regions defined for each Polygons
Mapping:
site.map <- ggmap(basemap)+ #basemap: terrain background (no labels) from google maps platform console, centered on avg lat/long of all sites
geom_polygon(data = gtnp_clip, aes(x = long, y = lat, group = group), color = "black", linetype = "dashed", fill = "darkgreen", alpha = 0.25, size = 0.6)+ #Adding clipped polygon of GTNP
geom_point(data = map_dat, aes(x = long, y = lat, fill = source, group = NULL), color = "black", size = 3, pch = 21)+ #add sites, color code by source
labs(color = "Hydrologic Source", fill = "Hydrologic Source")+ #updating legend title
scale_fill_manual(values = pres.pal, labels = c("Glacier", "Snow melt", "Subterranean ice"))+
theme_void()+ #removing lat/long from x and y axes
theme(legend.title = element_text(size = 9, face = "bold"), #adjusting legend title/text/plot title size, face, and position
legend.text = element_text(size = 9),
legend.position = c(0.84, 0.12), #Moving legend into the plot area
legend.background = element_rect(fill = "white", color = "black", size = 0.5), #adding box behind legend
legend.margin = margin(t=5, r=5, b=5, l = 5, unit = "pt") #Increasing margins on legend box
)+
ggrepel::geom_label_repel(data=map_dat, aes(x=long, y = lat, label=full_name, fill = source), size =3.25, position = position_dodge(width = 1), fontface = "bold", show.legend = FALSE)+ #adding labels for each stream, color-coded by source
geom_text_npc(aes(npcx = 0.84, npcy = 0.65, label = "Grand Teton \n National Park"), size = 4.5, fontface = 4, color = "black", hjust = 0.5, alpha = 0.25)+ #add label for Grand Teton NP
geom_text_npc(aes(npcx = 0.2, npcy = 0.8, label = "Targhee \n National Forest"), size = 4.5, fontface = 4, color = "black", hjust = 0.5, alpha = 0.25) #add lable for Targhee NF
site.map
First, group by source and re-calculate max, mean, min, and degree days:
summer_source <- stream_summer %>%
group_by(source, year)%>% # group by source type and year
summarise(t_xbar = mean(summer_xbar),
t_max = mean(max_xbar),
t_min = mean(min_xbar), #calculate mean daily temp, mean daily min and max
t_range = mean(range_xbar),
dd2.5_xbar = mean(dd_2.5),
dd5_xbar = mean(dd_5),
dd10_xbar = mean(dd_10),#calculate mean no. days with max temp >2.5, 5, and 10 C.
n_sites = length(unique(site)))%>% #calculate no sites for each year
filter(!is.na(year), !is.na(source)
#, !n_sites < 2
)%>%
mutate(source = case_when(source == "glacier"~"Glacier", source == "snowmelt"~"Snow melt", source == "sub_ice"~"Subterranean ice"))
## `summarise()` has grouped output by 'source'. You can override using the
## `.groups` argument.
Visualization:
source.means <- ggplot(summer_source, aes(x = source, y = t_xbar, fill = source))+
geom_boxplot()+
theme_classic()+
scale_fill_manual(values = pres.pal, labels = c("Glacier", "Snowmelt", "Subterranean Ice"))+
theme(legend.position = "none")+
labs(x = "", y = "Mean daily summer (July-Aug) temp., C")
source.max <- ggplot(summer_source, aes(x = source, y = t_max, fill = source))+
geom_boxplot()+
theme_classic()+
scale_fill_manual(values = pres.pal, labels = c("Glacier", "Snowmelt", "Subterranean Ice"))+
theme(legend.position = "none")+
labs(x = "Source", y = "Mean daily max. summer (July-Aug) temp., C")
source.min <- ggplot(summer_source, aes(x = source, y = t_max, fill = source))+
geom_boxplot()+
theme_classic()+
scale_fill_manual(values = pres.pal, labels = c("Glacier", "Snowmelt", "Subterranean Ice"))+
theme(legend.position = "none")+
labs(x = "", y = "Mean daily min. summer (July-Aug) temp., C")
source.means|source.max|source.min
Snow melt streams look to be significantly warmer and have higher max and mins - no surprise there. Check significance of these differences:
max.aov <- aov(t_xbar~source, data = summer_source)
summary(max.aov)
## Df Sum Sq Mean Sq F value Pr(>F)
## source 2 86.33 43.17 44.81 4.08e-08 ***
## Residuals 20 19.26 0.96
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
TukeyHSD(max.aov)
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = t_xbar ~ source, data = summer_source)
##
## $source
## diff lwr upr p adj
## Snow melt-Glacier 4.07427956 2.832756 5.315803 0.0000002
## Subterranean ice-Glacier 0.01385116 -1.271248 1.298950 0.9995901
## Subterranean ice-Snow melt -4.06042840 -5.345527 -2.775330 0.0000003
Snowmelt significantly warmer than sub-ice or glacier, no difference between sub-ice and glacier fed.
First, visualization - re-organize data for easier plotting:
source_plot <- as.data.frame(summer_source) %>%
dplyr::select(source, year, t_xbar, t_max, t_min, t_range)%>%
pivot_longer(-c(source, year), names_to = "measure", values_to = "temp")
Plotting:
ylab <- expression(bold("Stream Temperature " (degree*C)))
source.summer <- ggplot(source_plot, aes(x = year, y = temp, color = measure))+
geom_point()+
geom_line()+
geom_smooth(se = F, method = "lm", linetype = "dashed", size = 0.5)+
theme_classic()+
facet_wrap(~source, scales = "free_x")+
scale_color_manual(values = c("darkred","darkblue", "darkorange","darkgreen" )
, labels = c("Max. temp.", "Min. temp.", "Temp. range", "Mean temp")
)+
labs(x = "Date", y = ylab, color = "Measurement Type")+
theme(legend.position = "bottom",
axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
axis.title = element_text(size = 11, face = "bold"),
axis.text = element_text(size = 8),
legend.text = element_text(size = 10),
legend.title = element_text(size = 11, face = "bold"),
strip.text = element_text(size = 10, face = "bold"))
source.summer
## `geom_smooth()` using formula 'y ~ x'
Test for significance:
Nest plot data by source and measurement type:
source_nested <- source_plot%>%
nest(data = -c(source, measure))
Modeling:
source.lms <- source_nested %>%
mutate(model = purrr::map(data, ~lm(temp~year, data = .)%>% #run a model of temp~year for each site (referenced by "."); store results in new column called "model"
tidy))%>% #tidy function to construct a summary table with model results (seperate tables for each site)
unnest(model) %>% #unnest model column to get rows for each site, columns for model metrics
filter(term == "year") #remove intercept term to just get info on year for each site
source.lms
## # A tibble: 12 × 8
## source measure data term estimate std.error statistic p.value
## <chr> <chr> <list> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 Glacier t_xbar <tibble> year 0.139 0.0394 3.53 0.0124
## 2 Glacier t_max <tibble> year 0.0935 0.0722 1.30 0.243
## 3 Glacier t_min <tibble> year 0.146 0.0255 5.71 0.00125
## 4 Glacier t_range <tibble> year -0.0522 0.0606 -0.860 0.423
## 5 Snow melt t_xbar <tibble> year -0.189 0.239 -0.793 0.458
## 6 Snow melt t_max <tibble> year -0.149 0.291 -0.512 0.627
## 7 Snow melt t_min <tibble> year -0.182 0.201 -0.908 0.399
## 8 Snow melt t_range <tibble> year 0.0334 0.156 0.214 0.838
## 9 Subterranean ice t_xbar <tibble> year -0.0920 0.0928 -0.992 0.367
## 10 Subterranean ice t_max <tibble> year -0.167 0.202 -0.825 0.447
## 11 Subterranean ice t_min <tibble> year -0.0428 0.0584 -0.733 0.497
## 12 Subterranean ice t_range <tibble> year -0.124 0.191 -0.651 0.544
Significant increases in min, mean, and max summer temp in glacier-fed streams, no change in snowmelt fed streams, non-significant increases in sub_ice streams.
Calculate degrees of warming per decade for each:
warming_amounts <- source.lms %>%
mutate(deg_decade = 10*estimate)%>%
dplyr::select(source, measure, deg_decade, p.value)%>%
rename(Source = source, Measure = measure, Degrees.Warming = deg_decade, P = p.value)%>%
mutate(Measure = case_when(Measure == "t_xbar"~"Mean", Measure == "t_min"~"Minimum", Measure == "t_max"~"Maximum", Measure == "t_range"~"Diel Temp Change"))
warming_amounts
## # A tibble: 12 × 4
## Source Measure Degrees.Warming P
## <chr> <chr> <dbl> <dbl>
## 1 Glacier Mean 1.39 0.0124
## 2 Glacier Maximum 0.935 0.243
## 3 Glacier Minimum 1.46 0.00125
## 4 Glacier Diel Temp Change -0.522 0.423
## 5 Snow melt Mean -1.89 0.458
## 6 Snow melt Maximum -1.49 0.627
## 7 Snow melt Minimum -1.82 0.399
## 8 Snow melt Diel Temp Change 0.334 0.838
## 9 Subterranean ice Mean -0.920 0.367
## 10 Subterranean ice Maximum -1.67 0.447
## 11 Subterranean ice Minimum -0.428 0.497
## 12 Subterranean ice Diel Temp Change -1.24 0.544
Convert to FlexTable and export:
tbl2<-flextable(warming_amounts)%>% #convert to flextable
theme_vanilla()%>% #update theme to "vanilla"
set_table_properties(align = "center", layout = "autofit")%>% #center align the table on the screen, autofit columns
align(part = "all", align = "left")%>%#left align all text
set_header_labels(Source = "Hydrologic Source",
Measure = "Temperature Type",
Degrees.Warming = "Change per decade, C",
P = "P-value")
tbl2
Hydrologic Source | Temperature Type | Change per decade, C | P-value |
|---|---|---|---|
Glacier | Mean | 1.3909639 | 0.012412153 |
Glacier | Maximum | 0.9354068 | 0.242685598 |
Glacier | Minimum | 1.4571149 | 0.001248718 |
Glacier | Diel Temp Change | -0.5217080 | 0.422607515 |
Snow melt | Mean | -1.8949222 | 0.458180471 |
Snow melt | Maximum | -1.4886369 | 0.626870630 |
Snow melt | Minimum | -1.8226984 | 0.399083832 |
Snow melt | Diel Temp Change | 0.3340615 | 0.837940423 |
Subterranean ice | Mean | -0.9203200 | 0.366862680 |
Subterranean ice | Maximum | -1.6688166 | 0.447046968 |
Subterranean ice | Minimum | -0.4278540 | 0.496674074 |
Subterranean ice | Diel Temp Change | -1.2409626 | 0.543971003 |
save_as_image(tbl2, path = "manuscript/tables/degrees_warming.png", res = 300) #save as .png @ 300dpi
## [1] "manuscript/tables/degrees_warming.png"
save_as_docx(tbl2, path ="manuscript/tables/degrees_warming.docx", align = "center") #save as .png @ 300dpi
Adding these P-values to the plot:
source_pvals <- merge(source.lms, summer_source)%>% #add model info to source dataframe
mutate(ypos = case_when(measure == "t_max" ~ 0.95, measure == "t_xbar" ~ 0.9, measure == "t_min" ~ 0.85, measure == "t_range"~0.8)) #create column for adjusting P-value label Y position
source.summer.pv<- source.summer+
geom_text_npc(data = source_pvals, aes(npcx = 0.1, npcy = ypos, color = measure, label = paste( "P = ",round(p.value, 3), sep = ""))) # add p value labels to plot
source.summer.pv
## `geom_smooth()` using formula 'y ~ x'
Save:
ggsave("Temperature/plots/source_summer.pdf", source.summer.pv, device = "pdf", units = "in", height = 5, width = 6.5, dpi = "retina")
## `geom_smooth()` using formula 'y ~ x'
Re-organize data for easier plotting and modeling:
sourcedd_plot <- summer_source %>%
dplyr::select(source, year, dd2.5_xbar, dd5_xbar, dd10_xbar)%>%
pivot_longer(-c(source, year), names_to = "dd_type", values_to = "no_days")
Test for significance:
Nest plot data by source and measurement type:
sourcedd_nested <- sourcedd_plot%>%
nest(data = -c(source, dd_type))
Modeling:
sourcedd.lms <- sourcedd_nested %>%
mutate(model = purrr::map(data, ~lm(no_days~year, data = .)%>% #run a model of temp~year for each site (referenced by "."); store results in new column called "model"
tidy))%>% #tidy function to construct a summary table with model results (seperate tables for each site)
unnest(model) %>% #unnest model column to get rows for each site, columns for model metrics
filter(term == "year") #remove intercept term to just get info on year for each site
sourcedd.lms
## # A tibble: 9 × 8
## # Groups: source [3]
## source dd_type data term estimate std.error statistic p.value
## <chr> <chr> <list> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 Glacier dd2.5_x… <tibble> year 2.53 0.988 2.56 0.0428
## 2 Glacier dd5_xbar <tibble> year 1.20 0.473 2.54 0.0441
## 3 Glacier dd10_xb… <tibble> year 0 0 NaN NaN
## 4 Snow melt dd2.5_x… <tibble> year 2.74 1.26 2.17 0.0731
## 5 Snow melt dd5_xbar <tibble> year 2.68 1.27 2.11 0.0791
## 6 Snow melt dd10_xb… <tibble> year 2.05 1.30 1.57 0.167
## 7 Subterranean ice dd2.5_x… <tibble> year -0.475 1.38 -0.344 0.745
## 8 Subterranean ice dd5_xbar <tibble> year 0.0976 0.948 0.103 0.922
## 9 Subterranean ice dd10_xb… <tibble> year 0.0205 0.128 0.161 0.879
No significant changes in degree days in any group.
Plotting with these P-values:
sourcedd_pvals <- merge(sourcedd.lms, sourcedd_plot)%>% #add model info to source dataframe
mutate(ypos = case_when(dd_type == "dd2.5_xbar" ~ 0.85, dd_type == "dd5_xbar" ~ 0.9, dd_type == "dd10_xbar" ~ 0.95)) #create column for adjusting P-value label Y position
sourcedd.summer<- ggplot(sourcedd_pvals, aes(x = year, y = no_days, color = dd_type))+
geom_point()+
geom_line()+
geom_smooth(se = F, method = "lm", linetype = "dashed", size = 0.5)+
theme_classic()+
facet_wrap(~source, scales = "free_x")+
scale_color_manual(values = c("darkred","orange", "darkorange"), labels = c(">10 C", ">2.5 C", ">5 C"))+
labs(x = "Year", y = "Num. degree days", color = "Degree day threshold")+
theme(legend.position = "bottom",
axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
axis.title = element_text(size = 11, face = "bold"),
axis.text = element_text(size = 8),
legend.text = element_text(size = 10),
legend.title = element_text(size = 11, face = "bold"),
strip.text = element_text(size = 10, face = "bold"))+
geom_text_npc(aes(npcx = 0.9, npcy = ypos, color = dd_type, label = paste( "P = ",round(p.value, 3), sep = ""))) # add p value labels to plot
sourcedd.summer
## `geom_smooth()` using formula 'y ~ x'
First, combine glacier and snowmelt streams and recalculate mean, min, and max summer temperatures:
lumped_source <- stream_summer %>%
mutate(source_lumped = ifelse(source == "sub_ice", "sub_ice", "other"))%>% #add new col with lumped sources
filter(!is.na(year), !is.na(source))%>% #remove rows with no year or source
group_by(source_lumped, year)%>% #group by lumped source and year
summarise(t_xbar = mean(summer_xbar), #calculate mean, min, and max summer temps
max_xbar = mean(max_xbar),
min_xbar = mean(min_xbar))%>%
pivot_longer(!c(source_lumped, year), names_to = "measure", values_to = "temp") #pivot longer for plotting + analysis
## `summarise()` has grouped output by 'source_lumped'. You can override using the
## `.groups` argument.
Nest data by lumped source and measurement type:
lumped_nested <- lumped_source%>%
nest(data = -c(source_lumped, measure))
lm for each source + year combo:
lumped.lms <- lumped_nested %>%
mutate(model = purrr::map(data, ~lm(temp~year, data = .)%>% #run a model of temp~year for each site (referenced by "."); store results in new column called "model"
tidy))%>% #tidy function to construct a summary table with model results (seperate tables for each site)
unnest(model) %>% #unnest model column to get rows for each site, columns for model metrics
filter(term == "year") #remove intercept term to just get info on year for each site
lumped.lms
## # A tibble: 6 × 8
## # Groups: source_lumped [2]
## source_lumped measure data term estimate std.error statistic p.value
## <chr> <chr> <list> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 other t_xbar <tibble> year -0.0262 0.127 -0.206 0.843
## 2 other max_xbar <tibble> year -0.0300 0.180 -0.167 0.873
## 3 other min_xbar <tibble> year -0.0173 0.103 -0.168 0.872
## 4 sub_ice t_xbar <tibble> year -0.0920 0.0928 -0.992 0.367
## 5 sub_ice max_xbar <tibble> year -0.167 0.202 -0.825 0.447
## 6 sub_ice min_xbar <tibble> year -0.0428 0.0584 -0.733 0.497
Add p-values to data and plot:
lumped_pvals <- merge(lumped_source, lumped.lms)%>% #add model info to source dataframe
mutate(ypos = case_when(measure == "t_xbar" ~ 0.9, measure == "max_xbar" ~ 0.95, measure == "min_xbar" ~ 0.85)) #create column for adjusting P-value label Y position
lumped.summer<- ggplot(lumped_pvals, aes(x = year, y = temp, color = measure))+
geom_point()+
geom_line()+
geom_smooth(se = F, method = "lm", linetype = "dashed", size = 0.5)+
theme_classic()+
facet_wrap(~source_lumped, scales = "free_x")+
scale_color_manual(values = c("darkred","darkblue", "darkgreen"), labels = c("Max. temp.", "Min. temp.", "Mean temp"))+
labs(x = "Year", y = "Stream temperature, C", color = "Measurement Type")+
theme(legend.position = "bottom",
axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
axis.title = element_text(size = 11, face = "bold"),
axis.text = element_text(size = 8),
legend.text = element_text(size = 10),
legend.title = element_text(size = 11, face = "bold"),
strip.text = element_text(size = 10, face = "bold"))+
geom_text_npc(aes(npcx = 0.1, npcy = ypos, color = measure, label = paste( "P = ",round(p.value, 3), sep = ""))) # add p value labels to plot
lumped.summer
## `geom_smooth()` using formula 'y ~ x'
ggsave("Temperature/plots/lumped_summer.pdf", lumped.summer, device = "pdf", units = "in", height = 4, width = 6.5, dpi = "retina")
## `geom_smooth()` using formula 'y ~ x'
First, calculate April 1 SWE and summer air temp from snotel data; add to mean temperature dataset:
#Calc means across both sites:
snotel_means <- snotel %>%
mutate(date1 = as.Date(date, format = "%m/%d/%y"), #recode date column as R date
airtemp_c = as.numeric(airtemp_c))%>% #recode airtemp as numeric
group_by(date1)%>% #group by date
summarise(swe_xbar = mean(swe_cm), #calculate mean SWE and airtemp across both stations
airtemp_xbar = mean(airtemp_c, na.rm = T))%>%
mutate(mo = month(date1), yr = year(date1)) #Add month and year cols
## Warning: There was 1 warning in `mutate()`.
## ℹ In argument: `airtemp_c = as.numeric(airtemp_c)`.
## Caused by warning:
## ! NAs introduced by coercion
#Calculate max SWE and mean summer airtemp for each year:
snotel_summary <- snotel_means %>%
group_by(yr)%>% #group by year
mutate(swe_max = max(swe_xbar),
swe_max_date = date1[which.max(swe_xbar)],
swe_zero_date = date1[which.min(swe_xbar)][1])%>% #extract max SWE and add to new col
ungroup()%>% #ungroup to get all columns back
filter(mo == 7 | mo == 8)%>% #extract summer months
group_by(yr)%>% #group by year
summarise(airtemp_summer = mean(airtemp_xbar), #calculate mean summer air temp
swe_max = unique(swe_max),
swe_max_date = ymd(unique(swe_max_date)),
swe_zero_date =ymd(unique(swe_zero_date)))%>% #keep max SWE
mutate(melt_days = as.numeric(swe_zero_date-swe_max_date))%>%
rename(year = yr)
#Merge with long-term temperature averages:
stream_snotel <- merge(summer_3plus, snotel_summary)
LMs - first, nest by site:
stream_snotel_nested <- stream_snotel %>%
nest(data = - site)
lm for each site:
swe.temp.lms <- stream_snotel_nested %>%
mutate(model = purrr::map(data, ~lm(summer_xbar~swe_max, data = .)%>% #run a model of temp~year for each site (referenced by "."); store results in new column called "model"
tidy))%>% #tidy function to construct a summary table with model results (seperate tables for each site)
unnest(model) %>% #unnest model column to get rows for each site, columns for model metrics
filter(term == "swe_max") #remove intercept term to just get info on year for each site
swe.temp.lms
## # A tibble: 12 × 7
## site data term estimate std.error statistic p.value
## <chr> <list> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 s_cascade <tibble [5 × 18]> swe_max -0.0185 0.0454 -0.407 0.711
## 2 n_teton <tibble [8 × 18]> swe_max -0.0455 0.0332 -1.37 0.219
## 3 mid_teton <tibble [8 × 18]> swe_max 0.000589 0.00313 0.188 0.857
## 4 windcave <tibble [5 × 18]> swe_max -0.00725 0.00771 -0.941 0.416
## 5 s_ak_basin <tibble [6 × 18]> swe_max -0.00762 0.0101 -0.753 0.494
## 6 paintbrush <tibble [6 × 18]> swe_max -0.0477 0.0400 -1.19 0.299
## 7 delta <tibble [6 × 18]> swe_max -0.00137 0.00313 -0.438 0.684
## 8 silver_run <tibble [5 × 18]> swe_max -0.0324 0.00963 -3.37 0.0435
## 9 s_teton <tibble [6 × 18]> swe_max -0.0186 0.0311 -0.598 0.582
## 10 grizzly <tibble [5 × 18]> swe_max -0.0395 0.0488 -0.810 0.477
## 11 skillet <tibble [5 × 18]> swe_max -0.0211 0.0174 -1.22 0.311
## 12 cloudveil <tibble [4 × 18]> swe_max -0.0174 0.0191 -0.913 0.457
Add P values and merge:
swe_summer_ps <- merge(swe.temp.lms, stream_snotel)
swe.summer.temp <- ggplot(swe_summer_ps, aes(x = swe_max, y = summer_xbar, color = source))+
geom_point()+
geom_line()+
geom_smooth(se = F, method = "lm", linetype = "dashed", size = 0.5)+
theme_classic()+
facet_wrap(~full_name, scales = "free")+
scale_color_manual(values = pres.pal, labels = c("Glacier", "Snowmelt", "Icy Seep"))+
geom_text_npc(aes(npcx = 0.9, npcy = 0.9, label = paste( "P = ",round(p.value, 2), sep = "")))+
labs(x = "Annual maxium SWE, cm", y = "Mean daily July-Aug. stream temp., C", color = "Stream Source")+
theme(legend.position = "bottom",
axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
axis.title = element_text(size = 11, face = "bold"),
axis.text = element_text(size = 8),
legend.text = element_text(size = 10),
legend.title = element_text(size = 11, face = "bold"),
strip.text = element_text(size = 10, face = "bold"))
swe.summer.temp
## `geom_smooth()` using formula 'y ~ x'
lm for each site:
air.temp.lms <- stream_snotel_nested %>%
mutate(model = purrr::map(data, ~lm(summer_xbar~airtemp_summer, data = .)%>% #run a model of temp~year for each site (referenced by "."); store results in new column called "model"
tidy))%>% #tidy function to construct a summary table with model results (seperate tables for each site)
unnest(model) %>% #unnest model column to get rows for each site, columns for model metrics
filter(term == "airtemp_summer") #remove intercept term to just get info on year for each site
air.temp.lms
## # A tibble: 12 × 7
## site data term estimate std.error statistic p.value
## <chr> <list> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 s_cascade <tibble [5 × 18]> airtemp_su… -0.683 0.270 -2.53 0.0852
## 2 n_teton <tibble [8 × 18]> airtemp_su… -0.880 0.815 -1.08 0.322
## 3 mid_teton <tibble [8 × 18]> airtemp_su… -0.0492 0.0708 -0.695 0.513
## 4 windcave <tibble [5 × 18]> airtemp_su… 0.00876 0.240 0.0365 0.973
## 5 s_ak_basin <tibble [6 × 18]> airtemp_su… -0.405 0.305 -1.33 0.255
## 6 paintbrush <tibble [6 × 18]> airtemp_su… 0.110 1.73 0.0637 0.952
## 7 delta <tibble [6 × 18]> airtemp_su… 0.112 0.105 1.07 0.346
## 8 silver_run <tibble [5 × 18]> airtemp_su… -0.929 1.59 -0.585 0.600
## 9 s_teton <tibble [6 × 18]> airtemp_su… -0.123 1.20 -0.103 0.923
## 10 grizzly <tibble [5 × 18]> airtemp_su… -0.339 1.42 -0.238 0.827
## 11 skillet <tibble [5 × 18]> airtemp_su… 0.584 0.456 1.28 0.290
## 12 cloudveil <tibble [4 × 18]> airtemp_su… -0.224 0.451 -0.497 0.668
Add P values and plot:
air_summer_ps <- merge(air.temp.lms, stream_snotel)
swe.summer.temp <- ggplot(air_summer_ps, aes(x = airtemp_summer, y = summer_xbar, color = source))+
geom_point()+
geom_line()+
geom_smooth(se = F, method = "lm", linetype = "dashed", size = 0.5)+
theme_classic()+
facet_wrap(~full_name, scales = "free")+
scale_color_manual(values = pres.pal, labels = c("Glacier", "Snowmelt", "Icy Seep"))+
geom_text_npc(aes(npcx = 0.9, npcy = 0.9, label = paste( "P = ",round(p.value, 2), sep = "")))+
labs(x = "Mean daily July-Aug. air temp, C", y = "Mean daily July-Aug. stream temp., C", color = "Stream Source")+
theme(legend.position = "bottom",
axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
axis.title = element_text(size = 11, face = "bold"),
axis.text = element_text(size = 8),
legend.text = element_text(size = 10),
legend.title = element_text(size = 11, face = "bold"),
strip.text = element_text(size = 10, face = "bold"))
swe.summer.temp
## `geom_smooth()` using formula 'y ~ x'
lm for each site:
melt.temp.lms <- stream_snotel_nested %>%
mutate(model = purrr::map(data, ~lm(summer_xbar~melt_days, data = .)%>% #run a model of temp~year for each site (referenced by "."); store results in new column called "model"
tidy))%>% #tidy function to construct a summary table with model results (seperate tables for each site)
unnest(model) %>% #unnest model column to get rows for each site, columns for model metrics
filter(term == "melt_days") #remove intercept term to just get info on year for each site
melt.temp.lms
## # A tibble: 12 × 7
## site data term estimate std.error statistic p.value
## <chr> <list> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 s_cascade <tibble [5 × 18]> melt_days 0.00638 0.0768 0.0831 0.939
## 2 n_teton <tibble [8 × 18]> melt_days -0.144 0.0537 -2.68 0.0364
## 3 mid_teton <tibble [8 × 18]> melt_days -0.00883 0.00550 -1.61 0.160
## 4 windcave <tibble [5 × 18]> melt_days -0.0379 0.0185 -2.04 0.134
## 5 s_ak_basin <tibble [6 × 18]> melt_days 0.0308 0.0127 2.43 0.0721
## 6 paintbrush <tibble [6 × 18]> melt_days -0.0532 0.0828 -0.642 0.556
## 7 delta <tibble [6 × 18]> melt_days 0.000493 0.00599 0.0824 0.938
## 8 silver_run <tibble [5 × 18]> melt_days -0.0251 0.0296 -0.849 0.458
## 9 s_teton <tibble [6 × 18]> melt_days -0.101 0.0332 -3.05 0.0381
## 10 grizzly <tibble [5 × 18]> melt_days -0.0988 0.0476 -2.08 0.129
## 11 skillet <tibble [5 × 18]> melt_days -0.0184 0.0273 -0.675 0.548
## 12 cloudveil <tibble [4 × 18]> melt_days -0.00711 0.0242 -0.293 0.797
Add P values and plot:
melt_summer_ps <- merge(melt.temp.lms, stream_snotel)
melt.summer.temp <- ggplot(melt_summer_ps, aes(x = melt_days, y = summer_xbar, color = source))+
geom_point()+
geom_line()+
geom_smooth(se = F, method = "lm", linetype = "dashed", size = 0.5)+
theme_classic()+
facet_wrap(~full_name, scales = "free")+
scale_color_manual(values = pres.pal, labels = c("Glacier", "Snowmelt", "Icy Seep"))+
geom_text_npc(aes(npcx = 0.9, npcy = 0.9, label = paste( "P = ",round(p.value, 2), sep = "")))+
labs(x = "Mean daily July-Aug. air temp, C", y = "Mean daily July-Aug. stream temp., C", color = "Stream Source")+
theme(legend.position = "bottom",
axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
axis.title = element_text(size = 11, face = "bold"),
axis.text = element_text(size = 8),
legend.text = element_text(size = 10),
legend.title = element_text(size = 11, face = "bold"),
strip.text = element_text(size = 10, face = "bold"))
melt.summer.temp
## `geom_smooth()` using formula 'y ~ x'
First, add SNOTEL summary info to mean summer temps grouped by source:
source_snotel <- merge(summer_source, snotel_summary)
Re-organize data for analysis and plotting:
source_swe <- source_snotel %>%
dplyr::select(source, swe_max, t_xbar, t_min, t_max, t_range)%>% #dplyr::select temperature, source, and swe cols
pivot_longer(!c(source, swe_max), names_to = "measure", values_to = "temp") #pivot longer for plotting + analysis
Nest data by source and measure:
sourceSwe_nested <- source_swe %>%
nest(data = -c(source, measure))
LMs for each source and measure:
source.swe.lms <- sourceSwe_nested %>%
mutate(model = purrr::map(data, ~lm(temp~swe_max, data = .)%>% #run a model of temp~year for each site (referenced by "."); store results in new column called "model"
tidy))%>% #tidy function to construct a summary table with model results (seperate tables for each site)
unnest(model) %>% #unnest model column to get rows for each site, columns for model metrics
filter(term == "swe_max") #remove intercept term to just get info on year for each site
source.swe.lms
## # A tibble: 12 × 8
## source measure data term estimate std.error statistic p.value
## <chr> <chr> <list> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 Glacier t_xbar <tibble> swe_m… -0.00443 0.00821 -0.540 0.609
## 2 Glacier t_min <tibble> swe_m… -0.00477 0.00763 -0.626 0.554
## 3 Glacier t_max <tibble> swe_m… -0.00108 0.00992 -0.109 0.917
## 4 Glacier t_range <tibble> swe_m… 0.00369 0.00767 0.481 0.647
## 5 Snow melt t_xbar <tibble> swe_m… -0.0430 0.0250 -1.72 0.136
## 6 Snow melt t_min <tibble> swe_m… -0.0434 0.0191 -2.28 0.0629
## 7 Snow melt t_max <tibble> swe_m… -0.0354 0.0331 -1.07 0.326
## 8 Snow melt t_range <tibble> swe_m… 0.00807 0.0188 0.429 0.683
## 9 Subterranean ice t_xbar <tibble> swe_m… 0.00273 0.0131 0.208 0.844
## 10 Subterranean ice t_min <tibble> swe_m… 0.00956 0.00674 1.42 0.215
## 11 Subterranean ice t_max <tibble> swe_m… -0.0101 0.0276 -0.366 0.729
## 12 Subterranean ice t_range <tibble> swe_m… -0.0197 0.0242 -0.812 0.454
Add p-values and plot:
sourceSwe_pval<-merge(source_swe, source.swe.lms)%>%
mutate(ypos = case_when(measure == "t_xbar" ~ 0.9, measure == "t_max" ~ 0.95, measure == "t_min" ~ 0.85, measure == "t_range"~0.8)) #create column for adjusting P-value label Y position
source.swe<- ggplot(sourceSwe_pval, aes(x = swe_max, y = temp, color = measure))+
geom_point()+
geom_line()+
geom_smooth(se = F, method = "lm", linetype = "dashed", size = 0.5)+
theme_classic()+
facet_wrap(~source, scales = "free_x")+
scale_color_manual(values = c("darkred","darkblue", "darkorange","darkgreen"), labels = c("Max. temp.", "Min. temp.", "Temp. range", "Mean temp"))+
labs(x = "Annual maximum SWE, cm", y = "Stream temperature, C", color = "Measurement Type")+
theme(legend.position = "bottom",
axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
axis.title = element_text(size = 11, face = "bold"),
axis.text = element_text(size = 8),
legend.text = element_text(size = 10),
legend.title = element_text(size = 11, face = "bold"),
strip.text = element_text(size = 10, face = "bold"))+
geom_text_npc(aes(npcx = 0.9, npcy = ypos, color = measure, label = paste( "P = ",round(p.value, 3), sep = ""))) # add p value labels to plot
source.swe
## `geom_smooth()` using formula 'y ~ x'
Save:
ggsave("Temperature/plots/source_swe.pdf", source.swe, device = "pdf", units = "in", height = 4, width = 6.5, dpi = "retina")
## `geom_smooth()` using formula 'y ~ x'
Re-organize data for analysis and plotting:
source_air <- source_snotel %>%
dplyr::select(source, airtemp_summer, t_xbar, t_min, t_max, t_range)%>% #dplyr::select temperature, source, and swe cols
pivot_longer(!c(source, airtemp_summer), names_to = "measure", values_to = "temp") #pivot longer for plotting + analysis
Nest data by source and measure:
sourceAir_nested <- source_air %>%
nest(data = -c(source, measure))
LMs for each source and measure:
source.air.lms <- sourceAir_nested %>%
mutate(model = purrr::map(data, ~lm(temp~airtemp_summer, data = .)%>% #run a model of temp~year for each site (referenced by "."); store results in new column called "model"
tidy))%>% #tidy function to construct a summary table with model results (seperate tables for each site)
unnest(model) %>% #unnest model column to get rows for each site, columns for model metrics
filter(term == "airtemp_summer") #remove intercept term to just get info on year for each site
source.air.lms
## # A tibble: 12 × 8
## source measure data term estimate std.error statistic p.value
## <chr> <chr> <list> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 Glacier t_xbar <tibble> airte… 0.208 0.178 1.17 0.287
## 2 Glacier t_min <tibble> airte… 0.240 0.156 1.53 0.176
## 3 Glacier t_max <tibble> airte… 0.103 0.229 0.451 0.668
## 4 Glacier t_range <tibble> airte… -0.137 0.174 -0.784 0.463
## 5 Snow melt t_xbar <tibble> airte… -1.11 0.554 -2.01 0.0917
## 6 Snow melt t_min <tibble> airte… -0.955 0.470 -2.03 0.0882
## 7 Snow melt t_max <tibble> airte… -1.19 0.694 -1.71 0.138
## 8 Snow melt t_range <tibble> airte… -0.232 0.437 -0.532 0.614
## 9 Subterranean ice t_xbar <tibble> airte… -0.525 0.169 -3.10 0.0270
## 10 Subterranean ice t_min <tibble> airte… -0.156 0.161 -0.970 0.377
## 11 Subterranean ice t_max <tibble> airte… -1.15 0.337 -3.42 0.0189
## 12 Subterranean ice t_range <tibble> airte… -0.995 0.350 -2.84 0.0363
Add p-values and plot:
sourceAir_pval<-merge(source_air, source.air.lms)%>%
mutate(ypos = case_when(measure == "t_xbar" ~ 0.9, measure == "t_max" ~ 0.95, measure == "t_min" ~ 0.85, measure == "t_range"~0.8)) #create column for adjusting P-value label Y position
source.air<- ggplot(sourceAir_pval, aes(x = airtemp_summer, y = temp, color = measure))+
geom_point()+
geom_line()+
geom_smooth(se = F, method = "lm", linetype = "dashed", size = 0.5)+
theme_classic()+
facet_wrap(~source, scales = "free_x")+
scale_color_manual(values = c("darkred","darkblue", "darkorange","darkgreen"), labels = c("Max. temp.", "Min. temp.", "Temp. range", "Mean temp"))+
labs(x = "Annual maximum SWE, cm", y = "Stream temperature, C", color = "Measurement Type")+
labs(x = "Mean July-Aug. air temp., C", y = "Stream temperature, C", color = "Measurement Type")+
theme(legend.position = "bottom",
axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
axis.title = element_text(size = 11, face = "bold"),
axis.text = element_text(size = 8),
legend.text = element_text(size = 10),
legend.title = element_text(size = 11, face = "bold"),
strip.text = element_text(size = 10, face = "bold"))+
geom_text_npc(aes(npcx = 0.1, npcy = ypos, color = measure, label = paste( "P = ",round(p.value, 3), sep = ""))) # add p value labels to plot
source.air
## `geom_smooth()` using formula 'y ~ x'
Save:
ggsave("Temperature/plots/source_air.pdf", source.air, device = "pdf", units = "in", height = 4, width = 6.5, dpi = "retina")
## `geom_smooth()` using formula 'y ~ x'
Re-organize data for analysis and plotting:
source_melt <- source_snotel %>%
dplyr::select(source, melt_days, t_xbar, t_min, t_max, t_range)%>% #dplyr::select temperature, source, and swe cols
pivot_longer(!c(source, melt_days), names_to = "measure", values_to = "temp") #pivot longer for plotting + analysis
Nest data by source and measure:
sourceMelt_nested <- source_melt %>%
nest(data = -c(source, measure))
LMs for each source and measure:
source.melt.lms <- sourceMelt_nested %>%
mutate(model = purrr::map(data, ~lm(temp~melt_days, data = .)%>% #run a model of temp~year for each site (referenced by "."); store results in new column called "model"
tidy))%>% #tidy function to construct a summary table with model results (seperate tables for each site)
unnest(model) %>% #unnest model column to get rows for each site, columns for model metrics
filter(term == "melt_days") #remove intercept term to just get info on year for each site
source.melt.lms
## # A tibble: 12 × 8
## source measure data term estimate std.error statistic p.value
## <chr> <chr> <list> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 Glacier t_xbar <tibble> melt_… -0.00818 0.0173 -0.473 0.653
## 2 Glacier t_min <tibble> melt_… -0.0123 0.0157 -0.786 0.462
## 3 Glacier t_max <tibble> melt_… -0.00128 0.0208 -0.0614 0.953
## 4 Glacier t_range <tibble> melt_… 0.0111 0.0157 0.703 0.508
## 5 Snow melt t_xbar <tibble> melt_… -0.0751 0.0561 -1.34 0.230
## 6 Snow melt t_min <tibble> melt_… -0.0560 0.0495 -1.13 0.301
## 7 Snow melt t_max <tibble> melt_… -0.0937 0.0652 -1.44 0.201
## 8 Snow melt t_range <tibble> melt_… -0.0377 0.0369 -1.02 0.347
## 9 Subterranean ice t_xbar <tibble> melt_… 0.0138 0.0252 0.546 0.609
## 10 Subterranean ice t_min <tibble> melt_… 0.00267 0.0157 0.170 0.871
## 11 Subterranean ice t_max <tibble> melt_… 0.0290 0.0536 0.540 0.612
## 12 Subterranean ice t_range <tibble> melt_… 0.0263 0.0494 0.532 0.618
Add p-values and plot:
sourceMelt_pval<-merge(source_melt, source.melt.lms)%>%
mutate(ypos = case_when(measure == "t_xbar" ~ 0.9, measure == "t_max" ~ 0.95, measure == "t_min" ~ 0.85, measure == "t_range"~0.8)) #create column for adjusting P-value label Y position
source.melt<- ggplot(sourceMelt_pval, aes(x = melt_days, y = temp, color = measure))+
geom_point()+
geom_line()+
geom_smooth(se = F, method = "lm", linetype = "dashed", size = 0.5)+
theme_classic()+
facet_wrap(~source, scales = "free_x")+
scale_color_manual(values = c("darkred","darkblue", "darkorange","darkgreen"), labels = c("Max. temp.", "Min. temp.", "Temp. range", "Mean temp"))+
labs(x = "Annual maximum SWE, cm", y = "Stream temperature, C", color = "Measurement Type")+
labs(x = "Days from peak SWE to zero SWE", y = "Stream temperature, C", color = "Measurement Type")+
theme(legend.position = "bottom",
axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
axis.title = element_text(size = 11, face = "bold"),
axis.text = element_text(size = 8),
legend.text = element_text(size = 10),
legend.title = element_text(size = 11, face = "bold"),
strip.text = element_text(size = 10, face = "bold"))+
geom_text_npc(aes(npcx = 0.1, npcy = ypos, color = measure, label = paste( "P = ",round(p.value, 3), sep = ""))) # add p value labels to plot
source.melt
## `geom_smooth()` using formula 'y ~ x'
Save:
ggsave("Temperature/plots/source_meltout.pdf", source.melt, device = "pdf", units = "in", height = 4, width = 6.5, dpi = "retina")
## `geom_smooth()` using formula 'y ~ x'
General concept: want to test a linear mixed model of stream temperature (min, max, mean, and range) with:
Current approach will be to run seperate models for each different source type.
Calculate day of year for max SWE and zero SWE:
stream_snotel1 <- stream_snotel%>%
mutate(swe_max_doy =yday(swe_max_date),
swe_zero_doy=yday(swe_zero_date))
Extract predictor variables:
predictors <- stream_snotel1%>%
dplyr::select(year, airtemp_summer, swe_max, melt_days, swe_max_doy, swe_zero_doy)%>%
distinct()
Calculate correlation matrix and melt for plotting:
cormat <-round(cor(predictors), 2)
cormat[upper.tri(cormat)]<- NA
cormat_melt<-reshape2::melt(cormat, na.rm = T)
cormat_melt
## Var1 Var2 value
## 1 year year 1.00
## 2 airtemp_summer year 0.75
## 3 swe_max year -0.18
## 4 melt_days year -0.32
## 5 swe_max_doy year 0.58
## 6 swe_zero_doy year 0.22
## 8 airtemp_summer airtemp_summer 1.00
## 9 swe_max airtemp_summer 0.20
## 10 melt_days airtemp_summer -0.28
## 11 swe_max_doy airtemp_summer 0.72
## 12 swe_zero_doy airtemp_summer 0.42
## 15 swe_max swe_max 1.00
## 16 melt_days swe_max 0.20
## 17 swe_max_doy swe_max 0.26
## 18 swe_zero_doy swe_max 0.55
## 22 melt_days melt_days 1.00
## 23 swe_max_doy melt_days -0.68
## 24 swe_zero_doy melt_days 0.59
## 29 swe_max_doy swe_max_doy 1.00
## 30 swe_zero_doy swe_max_doy 0.19
## 36 swe_zero_doy swe_zero_doy 1.00
Plot correlation matrix:
ggplot(cormat_melt, aes(x = Var1, y = Var2, fill = value))+
geom_tile(color = "white")+
scale_fill_gradient2(low = "blue", high = "red", mid = "white",
midpoint = 0, limit = c(-1,1), space = "Lab",
name="Pearson\nCorrelation") +
theme_minimal()+
theme(axis.text.x = element_text(angle = 45, vjust = 1,
size = 10, hjust = 1),
axis.text.y = element_text(
size = 10),
axis.title = element_blank())+
coord_fixed()+
geom_text(aes(label = value))
Using a threshold of 0.5, correlations are:
Problematic variables are swe max/zero DOY, year. For now, dropping year (more relevant info contained in summer air temperature) and Swe max/zero DOY (info from both is roughly contained in melt_days)
Confirming that remaining variables don’t have correlations greater than 0.5:
cormat_filter <- cormat_melt%>%
filter(!(Var1 == "year"|Var1=="swe_max_doy"|Var1=="swe_zero_doy"),
!(Var2 == "year"|Var2=="swe_max_doy"|Var2=="swe_zero_doy"))
ggplot(cormat_filter, aes(x = Var1, y = Var2, fill = value))+
geom_tile(color = "white")+
scale_fill_gradient2(low = "blue", high = "red", mid = "white",
midpoint = 0, limit = c(-1,1), space = "Lab",
name="Pearson\nCorrelation") +
theme_minimal()+
theme(axis.text.x = element_text(angle = 45, vjust = 1,
size = 10, hjust = 1),
axis.text.y = element_text(
size = 10),
axis.title = element_blank())+
coord_fixed()+
geom_text(aes(label = value))
Data prep for plotting & modeling:
lmm_pivot <- stream_snotel1 %>%
dplyr::select(max_xbar, min_xbar, range_xbar, summer_xbar, source, airtemp_summer, swe_max, melt_days, site, year)%>%
pivot_longer(!c(source, airtemp_summer, swe_max, melt_days, site, year), names_to = "measure", values_to = "temp_c")%>%
mutate(transform = log10(temp_c))
#for plotting model results:
lmm_plot <- lmm_pivot%>%
pivot_longer(!c(source, site, year, measure, temp_c, transform), names_to = "predictor_name", values_to = "predictor_value")
Check distribution of response variables:
#without log10 transform:
ggplot(lmm_pivot, aes(x = temp_c))+
geom_density(fill = "grey", color = "black")+
facet_wrap(~measure)+
theme_classic()
#with log10 transform:
ggplot(lmm_pivot, aes(x = transform))+
geom_density(fill = "grey", color = "black")+
facet_wrap(~measure)+
theme_classic()
All are left-skewed, but a log transform just makes them left-skewed. Skipping log transform for now.
Load packages:
library(lme4)
## Loading required package: Matrix
##
## Attaching package: 'Matrix'
## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpack
##
## Attaching package: 'lme4'
## The following object is masked from 'package:raster':
##
## getData
library(lmerTest)
##
## Attaching package: 'lmerTest'
## The following object is masked from 'package:lme4':
##
## lmer
## The following object is masked from 'package:stats':
##
## step
Nest data by source:
lmm_nested <- lmm_pivot%>%
nest(data = -c(source, measure))
Run LMM for each source and measure - form:
(t_max, t_min, t_range, t_ave) ~ airtemp_summer + swe_max + melt_days + random(site)
max_temps <- filter(lmm_nested, measure == "max_xbar")
max.lmm <- max_temps %>%
mutate(model = purrr::map(data, ~lmerTest::lmer(temp_c~melt_days+airtemp_summer+swe_max+(1|site), data = .)%>% #run a model of temp~year for each site (referenced by "."); store results in new column called "model"
broom.mixed::tidy(effects="fixed", conf.int = TRUE)))%>% #tidy function to construct a summary table with model results (seperate tables for each site)
unnest(model)%>% #unnest model column to get rows for each site, columns for model metrics
mutate(p.value = round(p.value, 4))%>%
filter(!term == "(Intercept)")%>%
dplyr::select(source, term, estimate, std.error, statistic, df, p.value)
## boundary (singular) fit: see help('isSingular')
max.lmm
## # A tibble: 9 × 7
## source term estimate std.error statistic df p.value
## <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 sub_ice melt_days -0.0138 0.0421 -0.328 12.0 0.748
## 2 sub_ice airtemp_summer -1.18 0.429 -2.74 12.0 0.0179
## 3 sub_ice swe_max 0.00178 0.0184 0.0968 12.0 0.924
## 4 snowmelt melt_days -0.120 0.0375 -3.19 22.2 0.0042
## 5 snowmelt airtemp_summer -1.70 0.620 -2.74 22.6 0.0117
## 6 snowmelt swe_max -0.0219 0.0199 -1.10 22.5 0.282
## 7 glacier melt_days -0.00648 0.0155 -0.417 16.0 0.682
## 8 glacier airtemp_summer -0.0621 0.226 -0.275 16.1 0.787
## 9 glacier swe_max 0.00214 0.00886 0.242 16.1 0.812
Results:
Subterranean ice:
Significant negative correlation between summer air temp and max stream temp (maybe more ice melting?)
Snowmelt:
Significant negative correlation between melt length and max temp; same with summer air temp and max temp (??)
Glacier:
No significant effects.
Visualization:
ggplot(filter(lmm_plot, measure == "max_xbar"), aes(x = predictor_value, y = temp_c, color = source, group = site))+
geom_point(size = 0.5)+
geom_line(size = 0.5)+
geom_smooth(aes(group = NULL), method = "lm", se = T, linetype = "dashed", color = "black", size = 0.5)+
theme_classic()+
scale_color_manual(values = pres.pal, labels = c("Glacier", "Snowmelt", "Icy Seep"))+
facet_grid(source~predictor_name, scales = "free")+
theme(legend.position = "none")+
labs(x = "Predictor Value", y = "Mean maximum July-Aug stream temp, C")
## `geom_smooth()` using formula 'y ~ x'
Run LMM for each source and measure - form:
(t_max, t_min, t_range, t_ave) ~ airtemp_summer + swe_max + melt_days + random(site)
min_temps <- filter(lmm_nested, measure == "min_xbar")
min.lmm <- min_temps %>%
mutate(model = purrr::map(data, ~lmerTest::lmer(temp_c~melt_days+airtemp_summer+swe_max+(1|site), data = .)%>% #run a model of temp~year for each site (referenced by "."); store results in new column called "model"
broom.mixed::tidy(effects="fixed", conf.int = TRUE)))%>% #tidy function to construct a summary table with model results (seperate tables for each site)
unnest(model)%>% #unnest model column to get rows for each site, columns for model metrics
mutate(p.value = round(p.value, 4))%>%
filter(!term == "(Intercept)")%>%
dplyr::select(source, term, estimate, std.error, statistic, df, p.value)
min.lmm
## # A tibble: 9 × 7
## source term estimate std.error statistic df p.value
## <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 sub_ice melt_days -0.000123 0.00849 -0.0144 10.0 0.989
## 2 sub_ice airtemp_summer -0.0525 0.0888 -0.591 10.0 0.568
## 3 sub_ice swe_max -0.00722 0.00376 -1.92 10.0 0.0838
## 4 snowmelt melt_days -0.0684 0.0128 -5.34 22.2 0
## 5 snowmelt airtemp_summer -0.750 0.212 -3.54 22.5 0.0018
## 6 snowmelt swe_max -0.0221 0.00680 -3.24 22.5 0.0037
## 7 glacier melt_days -0.00525 0.00642 -0.818 16.0 0.426
## 8 glacier airtemp_summer 0.0540 0.0933 0.579 16.1 0.571
## 9 glacier swe_max -0.00654 0.00366 -1.79 16.1 0.0927
Results:
Subterranean ice:
Significant negative correlation between maximum SWE and min summer stream temp (Maybe less ice melting bc of snow insulation?)
Snowmelt:
Significant negative correlation between melt length, summer air temp, and swe max and summer stream temp. Melt length makes sense, others don’t really…
Glacier:
No significant effects.
Visualization:
ggplot(filter(lmm_plot, measure == "min_xbar"), aes(x = predictor_value, y = temp_c, color = source, group = site))+
geom_point(size = 0.5)+
geom_line(size = 0.5)+
geom_smooth(aes(group = NULL), method = "lm", se = T, linetype = "dashed", color = "black", size = 0.5)+
theme_classic()+
scale_color_manual(values = pres.pal, labels = c("Glacier", "Snowmelt", "Icy Seep"))+
facet_grid(source~predictor_name, scales = "free")+
theme(legend.position = "none")+
labs(x = "Predictor Value", y = "Mean minimum July-Aug stream temp, C")
## `geom_smooth()` using formula 'y ~ x'
Run LMM for each source and measure - form:
(t_max, t_min, t_range, t_ave) ~ airtemp_summer + swe_max + melt_days + random(site)
mean_temps <- filter(lmm_nested, measure == "summer_xbar")
mean.lmm <- mean_temps %>%
mutate(model = purrr::map(data, ~lmerTest::lmer(temp_c~melt_days+airtemp_summer+swe_max+(1|site), data = .)%>% #run a model of temp~year for each site (referenced by "."); store results in new column called "model"
broom.mixed::tidy(effects="fixed", conf.int = TRUE)))%>% #tidy function to construct a summary table with model results (seperate tables for each site)
unnest(model)%>% #unnest model column to get rows for each site, columns for model metrics
mutate(p.value = round(p.value, 4))%>%
filter(!term == "(Intercept)")%>%
dplyr::select(source, term, estimate, std.error, statistic, df, p.value)
mean.lmm
## # A tibble: 9 × 7
## source term estimate std.error statistic df p.value
## <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 sub_ice melt_days -0.000173 0.0187 -0.00927 10.1 0.993
## 2 sub_ice airtemp_summer -0.416 0.195 -2.13 10.1 0.0586
## 3 sub_ice swe_max -0.00406 0.00828 -0.491 10.1 0.634
## 4 snowmelt melt_days -0.0950 0.0215 -4.43 22.1 0.0002
## 5 snowmelt airtemp_summer -1.19 0.349 -3.41 24.2 0.0023
## 6 snowmelt swe_max -0.0219 0.0112 -1.95 23.5 0.0635
## 7 glacier melt_days -0.00520 0.00863 -0.602 16.0 0.556
## 8 glacier airtemp_summer 0.0124 0.126 0.0989 16.1 0.922
## 9 glacier swe_max -0.00429 0.00492 -0.871 16.1 0.397
Results:
Subterranean ice:
Significant negative correlation between summer air temp and mean summer stream temp (Maybe less ice melting bc of snow insulation?)
Snowmelt:
Significant negative correlation between melt length, summer air temp, and swe max and summer stream temp. Melt length makes sense, others don’t really…
Glacier:
No significant effects.
Visualization:
ggplot(filter(lmm_plot, measure == "summer_xbar"), aes(x = predictor_value, y = temp_c, color = source, group = site))+
geom_point(size = 0.5)+
geom_line(size = 0.5)+
geom_smooth(aes(group = NULL), method = "lm", se = T, linetype = "dashed", color = "black", size = 0.5)+
theme_classic()+
scale_color_manual(values = pres.pal, labels = c("Glacier", "Snowmelt", "Icy Seep"))+
facet_grid(source~predictor_name, scales = "free")+
theme(legend.position = "none")+
labs(x = "Predictor Value", y = "Mean daily July-Aug stream temp, C")
## `geom_smooth()` using formula 'y ~ x'
Run LMM for each source and measure - form:
(t_max, t_min, t_range, t_ave) ~ airtemp_summer + swe_max + melt_days + random(site)
temp_ranges <- filter(lmm_nested, measure == "range_xbar")
range.lmm <- temp_ranges %>%
mutate(model = purrr::map(data, ~lmerTest::lmer(temp_c~melt_days+airtemp_summer+swe_max+(1|site), data = .)%>% #run a model of temp~year for each site (referenced by "."); store results in new column called "model"
broom.mixed::tidy(effects="fixed", conf.int = TRUE)))%>% #tidy function to construct a summary table with model results (seperate tables for each site)
unnest(model)%>% #unnest model column to get rows for each site, columns for model metrics
mutate(p.value = round(p.value, 4))%>%
filter(!term == "(Intercept)")%>%
dplyr::select(source, term, estimate, std.error, statistic, df, p.value)
range.lmm
## # A tibble: 9 × 7
## source term estimate std.error statistic df p.value
## <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 sub_ice melt_days -0.00424 0.0432 -0.0982 10.2 0.924
## 2 sub_ice airtemp_summer -0.977 0.450 -2.17 10.5 0.0538
## 3 sub_ice swe_max 0.00448 0.0191 0.235 10.4 0.819
## 4 snowmelt melt_days -0.0508 0.0270 -1.89 22.1 0.0726
## 5 snowmelt airtemp_summer -0.977 0.447 -2.19 22.2 0.0397
## 6 snowmelt swe_max -0.00179 0.0143 -0.125 22.2 0.902
## 7 glacier melt_days -0.00113 0.0130 -0.0874 16.0 0.931
## 8 glacier airtemp_summer -0.114 0.189 -0.604 16.0 0.554
## 9 glacier swe_max 0.00861 0.00739 1.16 16.0 0.261
Results:
Subterranean ice:
Marginally significant negative correlation between summer air temp and temperature variability (warmer summers = less variability)
Snowmelt:
Significant negative correlation of melt length and summer air temp with temperature variability (longer melt and warmer summers = less variability). Seems a bit contradictory.
Glacier:
No significant effects.
Visualization:
ggplot(filter(lmm_plot, measure == "range_xbar"), aes(x = predictor_value, y = temp_c, color = source, group = site))+
geom_point(size = 0.5)+
geom_line(size = 0.5)+
geom_smooth(aes(group = NULL), method = "lm", se = T, linetype = "dashed", color = "black", size = 0.5)+
theme_classic()+
scale_color_manual(values = pres.pal, labels = c("Glacier", "Snowmelt", "Icy Seep"))+
facet_grid(source~predictor_name, scales = "free")+
theme(legend.position = "none")+
labs(x = "Predictor Value", y = "Mean July-Aug daily stream temp. change, C")
## `geom_smooth()` using formula 'y ~ x'
First, calculate richness, mean density, and total abundance for each site for each year; remove sites with only 1 year of data:
oneYr <- density_source%>%
group_by(site)%>% #group by site
summarise(nYrs = length(unique(Year)))%>% #count no. years
filter(nYrs <3) #extract sites w/ only 1-2 yrs data
invert_means <- density_source %>%
group_by(stream_code, full_name, site, Year)%>% #group by stream and year
summarise(richness = length(unique(taxa_lwr)), #calculate yearly richness, density, and total abundance
density_xbar = mean(density_xbar),
tot_abund = sum(abundance_total),
source = unique(source))%>%
filter(!site%in%oneYr$site) #drop sites with 1-2yrs of data
## `summarise()` has grouped output by 'stream_code', 'full_name', 'site'. You can
## override using the `.groups` argument.
First, nest data by site:
invert_nested <- invert_means %>%
nest(data = -site)
LM’s of Richness ~ year for each site:
richness.yr.lm <- invert_nested %>%
mutate(model = purrr::map(data, ~lm(richness~Year, data = .)%>% #run a model of temp~year for each site (referenced by "."); store results in new column called "model"
tidy))%>% #tidy function to construct a summary table with model results (seperate tables for each site)
unnest(model) %>% #unnest model column to get rows for each site, columns for model metrics
filter(term == "Year") #remove intercept term to just get info on year for each site
richness.yr.lm
## # A tibble: 12 × 7
## # Groups: site [12]
## site data term estimate std.error statistic p.value
## <chr> <list> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 s_ak_basin <gropd_df [7 × 7]> Year -0.893 0.724 -1.23 0.272
## 2 cloudveil <gropd_df [4 × 7]> Year 2.7 0.794 3.40 0.0766
## 3 delta <gropd_df [6 × 7]> Year -0.714 0.654 -1.09 0.336
## 4 grizzly <gropd_df [6 × 7]> Year 0.457 0.556 0.822 0.458
## 5 gusher <gropd_df [3 × 7]> Year 2.50 2.60 0.962 0.512
## 6 mid_teton <gropd_df [8 × 7]> Year -0.429 0.397 -1.08 0.322
## 7 n_teton <gropd_df [8 × 7]> Year -0.774 0.794 -0.975 0.367
## 8 paintbrush <gropd_df [6 × 7]> Year -0.286 0.379 -0.753 0.493
## 9 s_cascade <gropd_df [7 × 7]> Year -2.28 0.756 -3.02 0.0294
## 10 skillet <gropd_df [6 × 7]> Year 0.343 0.845 0.406 0.706
## 11 s_teton <gropd_df [8 × 7]> Year -0.774 0.717 -1.08 0.322
## 12 windcave <gropd_df [8 × 7]> Year -0.464 0.314 -1.48 0.190
Add p-values and plot:
richness_withps <- merge(invert_means, richness.yr.lm)
pres.pal <- c("#C26ED6", "#F89225", "#76D96F")
richness.yr <- ggplot(richness_withps, aes(x = Year, y = richness, color = source))+
geom_point()+
geom_line()+
geom_smooth(se = F, method = "lm", linetype = "dashed", size = 0.5)+
theme_classic()+
facet_wrap(~full_name, scales = "free")+
scale_color_manual(values = pres.pal, labels = c("Glacier", "Snowmelt", "Icy Seep"))+
geom_text_npc(aes(npcx = 0.9, npcy = 0.9, label = paste( "P = ",round(p.value, 2), sep = "")))+
labs(x = "Date", y = "Taxonomic richness", color = "Stream Source")+
theme(legend.position = "bottom",
axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
axis.title = element_text(size = 11, face = "bold"),
axis.text = element_text(size = 8),
legend.text = element_text(size = 10),
legend.title = element_text(size = 11, face = "bold"),
strip.text = element_text(size = 10, face = "bold"))
richness.yr
## `geom_smooth()` using formula 'y ~ x'
Significant decline in richness at S Cascade, non-significant declines at most other sites
Save:
ggsave("invert_data/plots/richness_year.pdf", richness.yr, device = "pdf", dpi = "retina", units = "in", height = 5, width = 6.5)
## `geom_smooth()` using formula 'y ~ x'
LM’s of Density ~ year for each site:
dens.yr.lm <- invert_nested %>%
mutate(model = purrr::map(data, ~lm(density_xbar~Year, data = .)%>% #run a model of temp~year for each site (referenced by "."); store results in new column called "model"
tidy))%>% #tidy function to construct a summary table with model results (seperate tables for each site)
unnest(model) %>% #unnest model column to get rows for each site, columns for model metrics
filter(term == "Year") #remove intercept term to just get info on year for each site
dens.yr.lm
## # A tibble: 12 × 7
## # Groups: site [12]
## site data term estimate std.error statistic p.value
## <chr> <list> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 s_ak_basin <gropd_df [7 × 7]> Year 17.3 8.99 1.92 0.113
## 2 cloudveil <gropd_df [4 × 7]> Year 64.2 7.57 8.49 0.0136
## 3 delta <gropd_df [6 × 7]> Year 17.3 25.8 0.670 0.540
## 4 grizzly <gropd_df [6 × 7]> Year -13.6 31.5 -0.430 0.689
## 5 gusher <gropd_df [3 × 7]> Year 28.2 16.7 1.69 0.340
## 6 mid_teton <gropd_df [8 × 7]> Year 23.1 31.6 0.730 0.493
## 7 n_teton <gropd_df [8 × 7]> Year -7.91 13.9 -0.571 0.589
## 8 paintbrush <gropd_df [6 × 7]> Year -55.1 116. -0.473 0.661
## 9 s_cascade <gropd_df [7 × 7]> Year 8.02 13.7 0.585 0.584
## 10 skillet <gropd_df [6 × 7]> Year 25.8 20.4 1.27 0.273
## 11 s_teton <gropd_df [8 × 7]> Year 2.79 1.61 1.73 0.134
## 12 windcave <gropd_df [8 × 7]> Year 16.8 16.1 1.04 0.337
Add p-values and plot:
density_withps <- merge(invert_means, dens.yr.lm)
pres.pal <- c("#C26ED6", "#F89225", "#76D96F")
dens.yr <- ggplot(density_withps, aes(x = Year, y = density_xbar, color = source))+
geom_point()+
geom_line()+
geom_smooth(se = F, method = "lm", linetype = "dashed", size = 0.5)+
theme_classic()+
facet_wrap(~full_name, scales = "free")+
scale_color_manual(values = pres.pal, labels = c("Glacier", "Snowmelt", "Icy Seep"))+
geom_text_npc(aes(npcx = 0.9, npcy = 0.9, label = paste( "P = ",round(p.value, 2), sep = "")))+
labs(x = "Date", y = "Mean density, individuals/m2", color = "Stream Source")+
theme(legend.position = "bottom",
axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
axis.title = element_text(size = 11, face = "bold"),
axis.text = element_text(size = 8),
legend.text = element_text(size = 10),
legend.title = element_text(size = 11, face = "bold"),
strip.text = element_text(size = 10, face = "bold"))
dens.yr
## `geom_smooth()` using formula 'y ~ x'
No significant trends in average density at any site
Save:
ggsave("invert_data/plots/density_year.pdf", dens.yr, device = "pdf", dpi = "retina", units = "in", height = 5, width = 6.5)
## `geom_smooth()` using formula 'y ~ x'
LM’s of Density ~ year for each site:
abund.yr.lm <- invert_nested %>%
mutate(model = purrr::map(data, ~lm(tot_abund~Year, data = .)%>% #run a model of temp~year for each site (referenced by "."); store results in new column called "model"
tidy))%>% #tidy function to construct a summary table with model results (seperate tables for each site)
unnest(model) %>% #unnest model column to get rows for each site, columns for model metrics
filter(term == "Year") #remove intercept term to just get info on year for each site
abund.yr.lm
## # A tibble: 12 × 7
## # Groups: site [12]
## site data term estimate std.error statistic p.value
## <chr> <list> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 s_ak_basin <gropd_df [7 × 7]> Year 120. 111. 1.09 0.327
## 2 cloudveil <gropd_df [4 × 7]> Year 528. 128. 4.14 0.0536
## 3 delta <gropd_df [6 × 7]> Year -30.9 255. -0.121 0.909
## 4 grizzly <gropd_df [6 × 7]> Year -10.8 168. -0.0641 0.952
## 5 gusher <gropd_df [3 × 7]> Year 469. 390. 1.20 0.442
## 6 mid_teton <gropd_df [8 × 7]> Year 17.5 141. 0.124 0.906
## 7 n_teton <gropd_df [8 × 7]> Year -22.0 81.4 -0.270 0.796
## 8 paintbrush <gropd_df [6 × 7]> Year -276. 202. -1.36 0.244
## 9 s_cascade <gropd_df [7 × 7]> Year -129. 63.6 -2.03 0.0977
## 10 skillet <gropd_df [6 × 7]> Year 233. 185. 1.26 0.277
## 11 s_teton <gropd_df [8 × 7]> Year -98.4 75.2 -1.31 0.239
## 12 windcave <gropd_df [8 × 7]> Year 20.7 90.6 0.228 0.827
Add p-values and plot:
abund_withps <- merge(invert_means, abund.yr.lm)
pres.pal <- c("#C26ED6", "#F89225", "#76D96F")
abund.yr <- ggplot(density_withps, aes(x = Year, y = tot_abund, color = source))+
geom_point()+
geom_line()+
geom_smooth(se = F, method = "lm", linetype = "dashed", size = 0.5)+
theme_classic()+
facet_wrap(~full_name, scales = "free")+
scale_color_manual(values = pres.pal, labels = c("Glacier", "Snowmelt", "Icy Seep"))+
geom_text_npc(aes(npcx = 0.9, npcy = 0.9, label = paste( "P = ",round(p.value, 2), sep = "")))+
labs(x = "Date", y = "Total abundance", color = "Stream Source")+
theme(legend.position = "bottom",
axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
axis.title = element_text(size = 11, face = "bold"),
axis.text = element_text(size = 8),
legend.text = element_text(size = 10),
legend.title = element_text(size = 11, face = "bold"),
strip.text = element_text(size = 10, face = "bold"))
abund.yr
## `geom_smooth()` using formula 'y ~ x'
No significant trends in total abundance at any site.
Save:
ggsave("invert_data/plots/abundance_year.pdf", abund.yr, device = "pdf", dpi = "retina", units = "in", height = 5, width = 6.5)
## `geom_smooth()` using formula 'y ~ x'
First, group by source and calculate mean richness, density, and abundance:
source_inverts <- invert_means %>%
group_by(source, Year)%>%
summarise(richness_xbar = mean(richness),
density_xbar = mean(density_xbar),
abund_xbar = mean(tot_abund))
## `summarise()` has grouped output by 'source'. You can override using the
## `.groups` argument.
First, nest data by source:
source_inverts_nested <- source_inverts %>%
nest(data = -source)
LMs of richness ~ year for each source:
ricness.source.lm <- source_inverts_nested %>%
mutate(model = purrr::map(data, ~lm(richness_xbar~Year, data = .)%>% #run a model of temp~year for each site (referenced by "."); store results in new column called "model"
tidy))%>% #tidy function to construct a summary table with model results (seperate tables for each site)
unnest(model) %>% #unnest model column to get rows for each site, columns for model metrics
filter(term == "Year") #remove intercept term to just get info on year for each site
ricness.source.lm
## # A tibble: 3 × 7
## # Groups: source [3]
## source data term estimate std.error statistic p.value
## <chr> <list> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 glacier <tibble [8 × 4]> Year -0.00794 0.303 -0.0262 0.980
## 2 snowmelt <tibble [8 × 4]> Year -1.10 0.751 -1.46 0.194
## 3 sub_ice <tibble [8 × 4]> Year -0.819 0.474 -1.73 0.134
Add p-values and plot:
richnessSource_pv <- merge(source_inverts, ricness.source.lm)
pres.pal <- c("#C26ED6", "#F89225", "#76D96F")
richness.source <- ggplot(richnessSource_pv, aes(x = Year, y = richness_xbar, color = source))+
geom_point()+
geom_line()+
geom_smooth(se = F, method = "lm", linetype = "dashed", size = 0.5)+
theme_classic()+
facet_wrap(~source)+
scale_color_manual(values = pres.pal, labels = c("Glacier", "Snowmelt", "Icy Seep"))+
geom_text_npc(aes(npcx = 0.9, npcy = 0.95, label = paste( "P = ",round(p.value, 2), sep = "")))+
labs(x = "Year", y = "Mean taxonomic richness", color = "Stream Source")+
theme(legend.position = "none",
axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
axis.title = element_text(size = 11, face = "bold"),
axis.text = element_text(size = 8),
legend.text = element_text(size = 10),
legend.title = element_text(size = 11, face = "bold"),
strip.text = element_text(size = 10, face = "bold"))
richness.source
## `geom_smooth()` using formula 'y ~ x'
Significant decline in richness at subterranean ice sites
Save:
ggsave("invert_data/plots/richness_source.pdf", richness.source, device = "pdf", dpi = "retina", units = "in", height = 5, width = 6.5)
## `geom_smooth()` using formula 'y ~ x'
LMs of density ~ year for each source:
dens.source.lm <- source_inverts_nested %>%
mutate(model = purrr::map(data, ~lm(density_xbar~Year, data = .)%>% #run a model of temp~year for each site (referenced by "."); store results in new column called "model"
tidy))%>% #tidy function to construct a summary table with model results (seperate tables for each site)
unnest(model) %>% #unnest model column to get rows for each site, columns for model metrics
filter(term == "Year") #remove intercept term to just get info on year for each site
dens.source.lm
## # A tibble: 3 × 7
## # Groups: source [3]
## source data term estimate std.error statistic p.value
## <chr> <list> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 glacier <tibble [8 × 4]> Year -8.24 20.2 -0.408 0.698
## 2 snowmelt <tibble [8 × 4]> Year -7.58 23.2 -0.327 0.755
## 3 sub_ice <tibble [8 × 4]> Year 11.2 8.50 1.32 0.235
Add p-values and plot:
densSource_pv <- merge(source_inverts, dens.source.lm)
pres.pal <- c("#C26ED6", "#F89225", "#76D96F")
dens.source <- ggplot(densSource_pv, aes(x = Year, y = density_xbar, color = source))+
geom_point()+
geom_line()+
geom_smooth(se = F, method = "lm", linetype = "dashed", size = 0.5)+
theme_classic()+
facet_wrap(~source)+
scale_color_manual(values = pres.pal, labels = c("Glacier", "Snowmelt", "Icy Seep"))+
geom_text_npc(aes(npcx = 0.9, npcy = 0.95, label = paste( "P = ",round(p.value, 2), sep = "")))+
labs(x = "Year", y = "Mean density, individuals/m2", color = "Stream Source")+
theme(legend.position = "none",
axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
axis.title = element_text(size = 11, face = "bold"),
axis.text = element_text(size = 8),
legend.text = element_text(size = 10),
legend.title = element_text(size = 11, face = "bold"),
strip.text = element_text(size = 10, face = "bold"))
dens.source
## `geom_smooth()` using formula 'y ~ x'
No significant trends in average density in any source.
Save:
ggsave("invert_data/plots/density_source.pdf", dens.source, device = "pdf", dpi = "retina", units = "in", height = 5, width = 6.5)
## `geom_smooth()` using formula 'y ~ x'
LMs of abundance ~ year for each source:
abund.source.lm <- source_inverts_nested %>%
mutate(model = purrr::map(data, ~lm(abund_xbar~Year, data = .)%>% #run a model of temp~year for each site (referenced by "."); store results in new column called "model"
tidy))%>% #tidy function to construct a summary table with model results (seperate tables for each site)
unnest(model) %>% #unnest model column to get rows for each site, columns for model metrics
filter(term == "Year") #remove intercept term to just get info on year for each site
abund.source.lm
## # A tibble: 3 × 7
## # Groups: source [3]
## source data term estimate std.error statistic p.value
## <chr> <list> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 glacier <tibble [8 × 4]> Year -40.6 124. -0.327 0.755
## 2 snowmelt <tibble [8 × 4]> Year -104. 63.4 -1.64 0.152
## 3 sub_ice <tibble [8 × 4]> Year 8.48 74.8 0.113 0.913
Add p-values and plot:
abundSource_pv <- merge(source_inverts, abund.source.lm)
pres.pal <- c("#C26ED6", "#F89225", "#76D96F")
abund.source <- ggplot(abundSource_pv, aes(x = Year, y = abund_xbar, color = source))+
geom_point()+
geom_line()+
geom_smooth(se = F, method = "lm", linetype = "dashed", size = 0.5)+
theme_classic()+
facet_wrap(~source)+
scale_color_manual(values = pres.pal, labels = c("Glacier", "Snowmelt", "Icy Seep"))+
geom_text_npc(aes(npcx = 0.9, npcy = 0.95, label = paste( "P = ",round(p.value, 2), sep = "")))+
labs(x = "Year", y = "Mean density, individuals/m2", color = "Stream Source")+
theme(legend.position = "none",
axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
axis.title = element_text(size = 11, face = "bold"),
axis.text = element_text(size = 8),
legend.text = element_text(size = 10),
legend.title = element_text(size = 11, face = "bold"),
strip.text = element_text(size = 10, face = "bold"))
abund.source
## `geom_smooth()` using formula 'y ~ x'
No significant trends in average abundance in any source.
Save:
ggsave("invert_data/plots/abundance_source.pdf", abund.source, device = "pdf", dpi = "retina", units = "in", height = 5, width = 6.5)
## `geom_smooth()` using formula 'y ~ x'